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ABSTRACT 

Contemporary galactic star formation occurs predominantly within gravitationally 
unstable, cold, dense molecular gas within supersonic, turbulent, magnetized giant 
molecular clouds (GMCs). Significantly, because the chemical evolution timescale and 
the turbulent eddy-turnover timescale are comparable at typical GMC conditions, 
molecules evolve via inherently non-equilibrium chemistry which is strongly coupled 
to the dynamical evolution of the cloud. 

Current numerical simulation techniques, which include at most three decades 
in length scale, can just begin to bridge the divide between the global dynamical 
time of supersonic turbulent GMCs, and the thermal and chemical evolution within 
the thin post-shock cooling layers of their background turbulence. We address this 
GMC astrochemical scales problem using a solution methodology, which permits both 
complex three-dimensional turbulent dynamics as well as accurate treatment of non- 
equilibrium post-shock thermodynamics and chemistry. 

We present the current methodology in the context of the larger scope of physical 
processes important in understanding the chemical evolution of GMCs, including gas- 
phase chemistry, dust grains, surface chemistry, and turbulent heating. We present 
results of a new Lagrangian verification test for supersonic turbulence. We charac- 
terize the evolution of these species according to the dimensionless local post-shock 
Damkohler number, which quantifies the ratio of the dynamical time in the post-shock 
cooling flow to the chemical reaction time of a given species. 

Lastly, we discuss implications of this work to the selection of GMC molecular 
tracers, and the zeroing of chemical clocks of GMC cores. 

Key words: astrochemistry, hydrodynamics, molecular processes, turbulence, stars: 
formation, ISM:molecules. 



1 INTRODUCTION 

Giant molecular clouds (GMCs) are rich dynamical struc- 
tures resulting from the interplay of many complex pro- 
cesses, including supersonic magnetized turbulence, self- 
gravity, chemical evolution, star formation, jets and out- 
flows, and radiative transfer. Because star formation occurs 
exclusively within gravitationally unstable cold dense gas 



within GMCs (Shu et al. 1987 McKee & Ostriker 20071, 



they play a crucial role in connecting larger-scale galactic 
dynamics to protostars, protostellar disks, and planet for- 
mation. The goal of this paper is to develop novel simulation 
techniques which will help elucidate the astrochemical evo- 
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lution of realistic, turbulent GMCs, and in turn yield insight 
into both observations and theoretical models of GMCs. 

GMCs are also host to a wide range of complex 
molecules, which play vital roles in both the diagnostics and 
the dynamics of the cloud. Observers have detected a wide 
range of molecular species in GMCs, including simple di- 
atomic molecules like CO and CS, up to complex organic 
molecules like formaldehyde (H2CO), methanol (ClfaOH) 
and amino acetonitrile, potentially a precursor to the amino 
acid glycine. Tracer molecules like CO and NH3 play a vital 
role as observational diagnostics of the density, temperature, 
and velocity dispersion of the cloud. Additionally, the chem- 
istry of key ions, such as HCO'*' , regulate the overall ioniza- 
tion level of the cloud, and establish the degree to which the 
magnetic field is coupled to the gas. 

Much of the work done to understand the chemical evo- 
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lution of GMCs assumes the hydrodynamic background of 
the cloud is either fixed or smoothly varying. Crucially, su- 
personic turbulence significantly modifies this quiescent pic- 
ture of the chemical structure of GMCs through two key 
mechanisms. Firstly, shock heating generates high temper- 
atures in post-shock cooling flow regions, and promotes gas 
phase reaction pathways with high activation energies - 
particularly of neutral-neutral chemistry - which are sup- 



pressed at lower temperatures |Padoan et al. 20001. Sec- 
ondly, parcels of gas are compressed and rarefied by the 
turbulence, resulting in a log-normal distribution of density 



values ( Padoan fc Nordlund|1997 Vazquez-Semadenl|| 1994 



Federrath et al. 2008 Passot & Vazquez-Semadeni 1998 1 



This broad distribution of densities results in significant de- 
partures of local reaction rates from the mean over the entire 
cloud ( [Clover fc Mac Low|2007a|b I . 

The Atacama Large Millimeter Array (ALMA) 
promises unprecedented observations of molecular line tran- 
sitions in the millimeter-range with spatial resolution of .01 
arcseconds (or roughly 1 AU at the distance of Taurus) - ten 
times better resolution than either Very Large Array (VLA) 
or the Hubble Space Telescope (HST), and high enough to 
begin to peer into the detailed shock dynamics of CMCs 
( Herbst||2008 1. The challenge to modelers is to match such 
high-quality CMC astrochemical observations with next- 
generation theoretical and computational tools which build 
upon and extend existing techniques. 

To date, the chemical evolution of GMCs has been mod- 
eled in one of two classes of approximations. The first set of 
approximations begin with extensive chemical networks and 
incorporate simplified dynamics - typically either spatially 
homogeneous or spherically-symmetric ( Aikawa et al.|2005 



2008[ ). The second includes a pioneering body of work which, 
for the first time, coupled realistic three-dimensional hy- 
drodynamical or magnetohydrodynamical turbulent simula- 
tions, with simplified chemical networks and molecular line 
cooling ([Pavlovski et aL 2002 2 0'06l|2004[ [Smith et al.|2004 



Glover fc Mac Low[[2007a|bi |Clover et al.[[2010[ ). Both ap- 
proaches have clear merits, but also significant limitations. 
Notably, as we will discuss in more detail below, current 
numerical simulation techniques are limited in their ability 
to resolve the thin post-shock cooling layers in turbulent 
GMCs. Ultimately what is needed to fully understand the 
astrochemical evolution of realistic GMCs is a full multi- 
dimensional simulation, coupled to extensive chemical net- 
works, and an accurate treatment of the thermodynamics of 
the gcis in thin cooling layers. As we discuss below, this goal 
presents significant technical challenges, which we directly 
address in this paper. 

Modeling the post-shock cooling region while also fully 
simulating the global dynamics of a turbulent supersonic 
clump presents a significant challenge to any 3D simulation, 
even with the power of adaptive mesh refinement. The post- 
shock cooling region behind a J-shock is of order ~ 1 AU in 
thickness at n = 10''' cm"'' for a shock velocity Vs ~ 10 km/s, 
and behind a C-shock of order ~ 10 AU thickness under 
similar conditions ( jShuU fc Hollenbach|1978| . In contrast, a 
supersonic CMC clump of lO'^M© is of order 1 pc in size, 
resulting in a dynamic range of four to five decades in length 
scale between the largest and smallest scales required to si- 
multaneously model both supersonic turbulence and molec- 



ular shock chemistry. This dynamic range is roughly one to 
two decades greater than even the most highly-resolved sim- 
ulations of turbul ent GMCs ( [Kritsuk et al.|2006[[2007| [Fed- 
errath et al.|2010j >. This large dynamic range of timescales is 
depicted in figure [l] which shows both the global chemical 
evolution of a representative species in a turbulent CMC 
clump along a single parcel, as well as its local hydrody- 
namic and chemical evolution through the post-shock flow 
of a single J-shock. 

Signiflcantly, for a typical supersonic CMC clump simu- 
lation with 3-D RMS Mach number 3.5 on an Eulerian 512^ 
mesh, the timestep is roughly 500 yr, assuming a Courant 
number of 0.5, implying the entire post-shock flow is only 
captured within ~ 2 timesteps. Consequently, the thickness 
of shocks is unresolved in multidimensional large-scale CMC 
simulations. This is a signiflcant limitation to understand- 
ing the global astrochemical evolution of a CMC, since non- 
equilibrium cooling and chemistry become significant in the 
post-shock cooling region where molecular line cooling de- 
termines the post-shock temperature. 

We address the CMC astrochemical scales problem us- 
ing a novel solution methodology, which permits both com- 
plex three-dimensional turbulent dynamics as well as an ac- 
curate treatment of non-equilibrium post-shock chemistry. 
The key idea is to separate the dynamical and chemical evo- 
lution into two weakly-coupled problems, and treat the wide 
range of scales between the dynamical time on the scale of a 
CMC clump tdyn ~ L/w ~ 10® yr to the cooling timescale in 
a thin post-shock layer, fcooi ~ 10^ yr. This decoupling of the 
dynamics from the molecular cooling and chemical evolution 
in the post-shock fiow is a valid starting point for an approx- 
imation because dust grain and molecular line cooling are 
highly efficient in establishing the temperature of the cloud, 
and the shock jump conditions themselves are simply fixed 
by conservation. Indeed, detailed three-dimensional simu- 
lations of turbulent GMCs including chemical, kinematic, 
and thermodynamic evolution in a CMC have demonstrated 
that the equation of state is well-approximated as isother- 
mal everywhere outside of the immediate radiative post- 



shock cooling layers (Pavlovski et al. 20061. Consequently, 



the large-scale supersonic hydrodynamical evolution can be 
accurately treated using an isothermal approximation to the 
Euler equations of hydrodynamics, while the detailed molec- 
ular and atomic cooling and chemical evolution in thin post- 
shock flows and in regions of intense shear dissipation can 
be treated separately in post-processing. 

This paper is organized as follows. In ^ we present 
a comprehensive overview of the primary physical processes 
which govern the physical and chemical evolution of a CMC, 
including gas-phase chemistry, dust grains, surface chem- 
istry, and turbulent heating. In i|3[ we describe in detail 
the numerical methodology developed to treat the physical 
and chemical evolution within supersonic isothermal tur- 
bulent GMCs, including gas dynamics, driven turbulence, 
Lagrangian tracer particles, gas-phase chemistry, gas-phase 
cooling, and J- shock modeling. Section |3.1.2| presents the 
dimensionless parameterization selected for our hydrody- 
namics models. Section !|4| we present new verification tests 
developed to test the new aspects of the numerical methods 
here. In section ^ , we present preliminary results for the 
astrochemical evolution along a sample Lagrangian trajec- 
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tory as well as a ensemble of Lagrangian trajectories cal- 
culated with our new numerical methods. We calculate the 
joint PDF of density and temperature, and infer the mass 
filling fraction of the warm molecular phase generated by 
turbulent heating within our models. Lastly, in section i|6j 
we discuss our findings and conclude. 



2 PHYSICAL PROCESSES 



2.1 Gas-Phase Chemistry 



In a key paper, iHerbst & Klemperer ( 1973 1 demonstrated 



that ion-neutral chemistry could provide formation channels 
for a number of molecular species, such as OH and H2O. 
Specifically, cosmic rays are able to ionize molecular hydro- 
gen even deep in the interior of GMCs, and thereby initiate 
ion-neutral chemical reactions. The classical Langevin rate 
of these ion-neutral reactions is independent of velocity, and 
therefore able to proceed even at the low mean background 
temperatures within GMCs. Consequently, the ion-neutral 
gas-phase molecular chemistry timescale is set by cosmic 
ray injection in dense GMCs and is independent of density 
( |Herbst fc Klemperer|p73l [Bergin fc Tafalla|[2007l ) , with a 
characteristic chemical time scale tchem yr. 

In addition to ion-neutral chemistry, neutral-neutral 
chemistry is also now understood to play an important role 
in the chemical evolution of GMCs even at low tempera- 
tures. Direct measurements of some neutral-neutral reac- 
tion rate have demonstrated they they have relatively rapid 
reaction rates ( > 10~^° cm'^ molecule"^ s~^) even at low 



temperatures down to 10 K ( Smith|1997 1. In addition, other 
neutral-neutral reactions can be activated in warm regions 
with sufficient gas kinetic temperature to overcome activa- 
tion barriers. Our initial models include only a few dozen 
simple carbon, nitrogen, and oxygen- bearing species (in- 
cluding, among others, C, CH+, CO, HCO+, O, O2, OH, 
H2O, N, and NH3 - see table|3]for a complete list) for which 
the rate reactions are well-determined. 

A realistic treatment of both ion-neutral and neutral- 
neutral gas phase chemistry is at the heart of this paper. 
To further the state-of-the-art in astrochemical modeling, 
we have developed a new astrochemistry framework, AS- 
TROCHEM, which builds upon and extends existing astro- 
chemical codes and databases using a sophisticated auto- 
matic chemical network pre-processor, as discussed in more 
detail in our methodology section (p^. ASTROCHEM al- 
lows for the efficient and flexible solution of large chemical 
networks, as discussed below in our methodology section. 
The chemical reaction network presently included in AS- 
TROCHEM have been tabulated in appendix[Aj along with 
the mathematical expression of reaction rate coefficients of 
reactions and their references. 



2.2 Dust Grains and Surface Chemistry 

In addition to the gas phase molecular species, dust grains 
also play a crucial role in the chemistry and thermodynam- 
ics of GMCs. Because hydrogen has a high surface mobil- 
ity on dust grains, a key feature of dust grain chemistry 
is the hydrogenation of molecules, providing an important 



formation mechanism for complex organic molecules includ- 
ing methanol ( van Dishoeck fc Blake|1998 l. Although there 
is a significant body of both theoretical and experimental 
work on the surface chemistry at conditions appropriate for 
GMCs (e.g., recent review by |Herbst et al.|2005| among oth- 
ers), the complexity of the subject is significant. In partic- 
ular, the reaction rates depend directly on the surface mo- 
bility of hydrogen and other light molecules, which in turn 
depends on complex interfacial physics, including whether 
the surface molecules are strongly bound to the surface 
through chemisorption or weakly bound through physisorb- 
tion. There are also outstanding research issues with regard 
to the methodology used to evolve surface chemical species; 
while rate equation methods are still widely used, this ap- 
proach formally breaks down for trace species for which the 
mean abundance may be less than a single molecule per 
grain. 

The chemistry models present in this paper consist 
solely of gas-phase chemistry, and include the effect of the 
dust grains on the molecular hydrogen chemistry and gas 
heating/cooling only. This model will later be extended to 
incorporate a grain surface chemistry reaction network (e.g., 
Herbst et al, 2005) using a rate reaction formalism built 
upon ASTROCHEM. 



2.3 Turbulent Heating 

Turbulence has long been recognized as playing an impor- 
tant role in the dynamics of GMCs, at least to Larson's sem- 
inal paper on the internal velocity dispersion in molecular 
clouds ( Larson||l981 1. Turbulence decays on the order of a 
dynamical time ( Low|199^ I , thereby leading to a significant 
source of heating for GMCs. This idea was initially explored 
in the context of atomic HI clouds. Earlier work has dis- 
cussed the impact of subsonic turbulence on the gas kinetic 
temperature of atomic HI clouds (Falgarone fc Pugct 1995). 
Extensions of this work using a simplified model of intermit- 
tent turbulent heating demonstrated that localized hot re- 
gions could form within otherwise cold HI clouds, thereby ac- 
tivating neutral-neutral gas phase chemical reactions which 



were suppressed at the background temperature (Falgarone 
et al.|1995| ). 



More recently, simulations of molecular gas in turbu- 
lent GMCs have demonstrated that turbulent heating alone 
can largely account for the mean temperature within GMCs 
(Pan & Padoan|2009 ). Because of the intermittent nature of 
the turbulent heating, post-shock temperatures in excess of 



1000 K can be reached in simulation models ( Pavlovski et al 



2006 1. Strong shears may also lead to a large dissipation rate; 



the specific turbulent heating rate per unit mass is relatively 
uncorrelated with density ( Pan fc Padoan|[2b09p . However, 
the majority of the mass of the turbulent GMC model is 
concentrated within shocks. Extensions of the She-Leveque 
model of turbulent intermittency reveal the dimension of the 
most dissipative structures in supersonic turbulence which 
best fit both numerical and observational data to be corre- 
sponding to shocks, and not corresponding to vortices, as 
is the case with subsonic turbulence ( Boldyrev|2002 Pan fc 
|Padoan|2009[ ) . Recent work has shed more light on intermit- 
tency models and fractal dimensions of the most dissipative 
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structures of supersonic turbulence (Schmidt et al. 2009 
Konstandin et al.||2012[ ). 



A major goal of the current paper is to devise a consis- 
tent numerical methodology which allows us to treat both 
turbulent shock heating and shear dissipation in supersonic 
isothermal GMCs with high temporal and spatial resolu- 
tion. As we will describe, our numerical simulations capture 
the shocks on the mesh, allowing us to include the effects 
of shock heating on subgrid scales using post-processing of 
Lagrangian tracer particles. Furthermore, we also take into 
account dissipation of shear energy on resolved scales. We 
discuss the methodological procedures which allow us to in- 
corporate these effects more fully in section ^ 



2.4 Characteristic Length Scales and 
Dimensionless Numbers 

We begin by identifying the key length and time scales rele- 
vant to the supersonic chemical evolution of the GMC. Key 
physical insight into the physical and chemical evolution of 
a GMC can be obtained by examining the characteristic 
time scales of the GMC, as well as the dimensionless num- 
bers characterizing its properties. Whereas recent work has 
elucidated the dimensionless numbers important in passive 



scalar mixing in isothermal supersonic turbulence (Pan & 
Scannapieco||20l"0 1, we focus on the characterization of the 
chemical evolution of the GMC. 

The Reynolds number Re quantifies the ratio of inertial 
to viscous forces, and is defined on the length scale I as: 



Re^!^^7.49xl0«f-i 
1 pc 



f " 

V1.22 kms-i 



5.00 X 1014 cm2 s-i 



(1) 



where v is the turbulent RMS velocity, and v is the kine- 
matic viscosity due to neutral-neutral molecular collisions 
(defined in equation B4 in appendix [b|. The Reynolds num- 



ber within GMCs is enormously larger than critical value for 
the transition to fully-developed turbulence measured in ter- 
restrial experiments ( Pope|[2"000 1, and so we expect GMCs 
to be fully turbulent. A wide variety of observational diag- 
nostics, including power spectra measured in both the veloc- 
ity and density fields, are consistent with this expectation 



( Elmegreen fc Scalo|2004a|b[ [Mac Low fc Klessen||2004 i. 

The Kolmogorov scale rjK is the length scale where the 
turbulent kinetic energy is ultimately dissipated due to vis- 
cous forces. Mathematically, 



3 \ 4 



3.49 X 10"^AU I 



.5.00 X 10" cm2 



2.02 X 10-4 erg cm-^g-is" 



where e is the specific turbulent kinetic energy dissipation 
rate, and v is the kinematic viscosity. For simplicity, we have 
assumed that e is independent of the length scale; here it 
has been calculated at the integral scale using fiducial val- 
ues for a representative Mach 3.5 clump. Equivalently, the 
Kolmogorov scale is the length scale where the Reynolds Re 
number is one. The eddy turnover time at the Kolmogorov 



scale is the local Kolmogorov time scale tK, 



tK 



(a 



29.2 years x 



.5.00 X 1014 cm2 s-i 



5.89 X 10-4 erg cm-^g-is-i 



(3) 



Both rjK and tK are vastly smaller than the computational 
grid size Ax and timestep At, respectively, in 3-D simu- 
lations. Consequently, any direct numerical simulation ap- 
proach to turbulent mixing is computationally prohibitive, 
and some effective modelling is required. 

The Knudsen (Kn) is the ratio between molecular mean 
free path and a representative physical length. This deter- 
mines the validity of the continuum hypothesis. If the Knud- 
sen number is << 1 , then the mean- free path length is sig- 
nificantly less than the typical flow scale and the continuum 
approximation is valid. The mean-free path length of the ion 
and neutral species are, however, different, so their Knudsen 
numbers vary. The neutral-neutral Knudsen number is given 
by: 



KUnn = 1.37 X 10" 



6.07 X 10-15 cm2 



4 X 104 cm-3 

L 
1 pc 



(4) 



The Kolmogorov length scale sets the smallest length rel- 
evant to turbulent flows. For a Kolmogorov length scale of 
~ 4.6 X 10~^ AU, the neutral-neutral Knudsen number Knnn 
is ~ 6.1 X 10"'^. Consequently the fluid approximation re- 
mains valid throughout the turbulent flow, even far beneath 
the length scales currently probed by numerical simulations. 

Next, we compare the turbulent eddy turnover 
timescale with the timescale required for ion-neutral molec- 
ular chemistry. In reactive chemical flows, the Damkohler 
number Da is defined as the ratio of teddy (0 ^o tchom- For Da 
^ 1, chemical equilibrium is reached on an essentially frozen 
background flow. In contrast, for Da 5C 1 (the "well-mixed 
reactor" regime), a turbulent eddy can turn over before ion- 
neutral chemical equilibrium is established. Consequently, 
for Da ^ 1, turbulence plays an essential role in determining 
the flnal chemical products. We define the turbulent eddy 
turnover time timescale teddy as: 

1/2 



eddy 



I 



IL 



(5) 



where we have made use of Larson's law, v{l) — vo{l/ L)^^^ , 
which specifies the scaling of the turbulent RMS velocity Vr 



with the length scale I (Larson 19811. Here vo is the RMS 
velocity on a characteristic scale L; we adopt vq = 1 km/s 
at L = 1 pc as fiducial values for galactic GMCs. 

The ion-neutral gas-phase molecular chemistry 
(2) timescale is set by cosmic ray injection in dense GMCs 
and is independent of density ( ,Bergin fc Tafalla 20071 - 
tchem ~ 3 X 10^ yr. 

For these fiducial values of galactic GMCs, the 
Damkohler number established by turbulence transitions 
through Da = 1 at a characteristic scale of roughly 1 pc. 



Using Larson's mean density-size relation, and scaling to a 
fiducial mean density of 10^ cm~^ at 0.1 pc, we find that this 
characteristic length scale corresponds to a mean density of 
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10 cm~ , which is characteristic of the density of large- 
scale clumps in galactic GMCs. On smaller length scales, 
the chemical evolution timescale and the eddy turnover 
timescale are comparable, and molecules evolve via inher- 
ently non-equilibrium chemistry which is strongly coupled 
to the dynamical evolution of the cloud. Therefore on the 
clump scale or smaller, the timescale for turbulent mixing 
is comparable to the ion-neutral chemical timescale. Conse- 
quently, on the clump scale or smaller, dynamical processes 
play a essential role in determining GMC chemical abun- 
dances. 



3 METHODOLOGY 

In the present work, we model the turbulence within 
the GMC as purely hydrodynamic, or equivalently, that 
the clouds are magnetically supercritical. Consequently, all 
shocks are modeled as J-shocks. While there is indirect the- 
oretical evidence that GMCs are magnetically supercritical 
with respect to their mean magnetic field ( |McKee fc Ostriker] 
20071, there is an absence of direct measurements constrain- 



ing the mass-to-flux ratio on large scales. Future work will 
build upon the current framework to incorporate the mag- 
netic field and ambipolar diffusion, and treat C-shocks. 

In addition to the effects of shock heating, we also take 
into account the solenoidal heating due to dissipation within 
vortices and strong shears. In contrast to previous work 
which includes subgrid modeling of this effect (Falgarone 



& Puget 1995 Pan & Padoan 20091, we account for the 



solenoidal heating on resolved scales only. Even on these re- 
solved scales, the effect of solenoidal heating is significant, 
and accounts for two-thirds of the total turbulent heating in 



our models {[ 3.7 1. 

An extensive body of simulations has demonstrated 
that undriven, supersonic turbulence decays rapidly on a dy- 



namical timescale (Low et al. 1998 1. Furthermore, the statis- 



tical properties of observed giant molecular clouds cores have 
been demonstrated to be in better agreement with driven, 
steady-state turbulence than with undriven, decaying turbu- 
lence ( Offner et al.|2008 1 . Consequently, we model turbulent 
giant molecular clouds as being driven over a narrow band 
of wavenumbers. On smaller scales, an inertial range is es- 
tablished. Recent work show that forcing of the turbulence 
significantly affects the density and velocity statistics of su- 
personic turbulence in both the inertial and injection range 



( Schmidt et al.|2009[[Federrath et al.|2010[[Konstandin et al. 
|2012[ ). Since we post-process Lagrangian tracers by detecting 
shocks to get chemical evolution in our paper, any change 
in density scaling will change the way our chemical network 
evolves. 

Our solution methodology consists of several fundamen- 
tal ingredients. The first involves Eulerian hydrodynamics 



(i 3.1 1. The Eulerian hydrodynamical flow is driven on large 



scales to establish a turbulent power spectrum {[3.2). The 



Eulerian calculation also has embedded within it Lgrangian 



tracer particles {[ 3.3 1 which follow the flow and capture 
the Lagrangian state of the system for later detailed gas- 
phase chemistry (|3.4[), and radiative shock chemistry post- 



processing 3.6 1 



3.1 Gas Dynamics 

We solve the Euler equations of hydrodynamics with driving 
source source terms, 



dp 
dt 



+ V • ipv) = 0, 



d{pv) 
dt 

djpE) 
dt 

P = (7 - l)pe. 



-f- V ■ (pvv) = -VP -I- /, 
+ V ■[{pE + P)v\=v- f, 



(6) 

(7) 

(8) 
(9) 



Here p is the mass density, v is the velocity, P is the 
gas pressure, E is the total energy density (related to the 
internal energy per unit mass, e, hy E — e + ^v^). f is a 
force density representing the forcing of the fluid on large 



scales to sustain driven turbulence against decay (see [3.2). 

We model isothermal driven turbulence with an adia- 
batic equation of state with a ratio of speciflc heats 7 ~ 1. 
An isothermal gas can be thought of as an adiabatic gas 
in the limit that the number of internal degrees of freedom 
approaches inflnity, and can be approximately modeled by 
taking a large but finite number of degrees of freedom. This 
approximation is easily accommodated by many existing hy- 
drodynamics solvers, and has a domain of validity provided 
that the total heat energy added to the system over the 
duration of the simulation is much less than the initial in- 
ternal energy. However, the specific turbulent heating rate, 
which scales as v'^ /{v/L) ~ / L ~ M"^ , places a strict 
upper-bound constraint on the choice of 7 for supersonic 
turbulence. 

Consider a supersonic turbulence simulation modelled 
using an adiabatic equation of state in a periodic domain. 
The initial specific internal energy is fixed by the isothermal 
sound speed and the ratio of specific heats, and is simply 

fiint 



4o/(7 



Over each dynamical time, the decay 



of turbulent energy will deposit ~ 1/2M c^^^ of heat en- 
ergy per unit mass into the system. Physically, this energy 
is radiated from an isothermal system, but under the adia- 
batic approximation, the energy is simply deposited as heat. 
Consequently, in order for the gas to remain approximately 
isothermal over a total evolution of A'dyn dynamical times, 
we require 



7- K 



(10) 



In the driven simulations presented here, one begins 
with quiescent gas and must evolve for a minimum of 
-^dyn ~ 2 — 3 to reach a steady-state. Typically, in order 
to gather sufficient turbulent statistics in steady-state, one 
wishes to run for A'dyn ~ 10. We find that for modelling 
a typical turbulent GMC clump-sized region, with = 5, 
and Adyn ~ 10, we require 7 — 1 <C .008. In the simulations 
presented here, we utilize 7 = 1-1- 10~^, which produces 
excellent steady-state results. While higher values of 7 are 
found in literature (often as large as 7 = 1 -I- lO"'^), these 
will not be sufficient to maintain isothermality throughout 
the course of a highly-supersonic simulation, and will lead 
to non-steady turbulence statistics. 
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3.1.1 Self-Consistent Choice of Mean Temperature and 
Mach Number 

The Eulerian frame calculations are completed using isother- 
mal hydrodynamics, and can be characterized by a single di- 
mensionless number, the 3-D RMS turbulent Mach number. 
The thermal evolution of Lagrangian fluid elements must 
further incorporate an energy equation including both heat- 
ing and cooling, and requires an initial temperature. After 
passing through shocks, Lagrangian fluid elements return to 
the initial temperature. An additional self-consistency re- 
quirement is therefore that the mean temperature within 
the Eulerian volume is equal to the mean temperature along 
Lagrangian particles. This is equivalent to asserting that the 
overall system is in net thermal equilibrium, once all cool- 
ing and heating processes, including turbulent heating, are 
taken into account. 

In order to set the mean temperature self-consistently 
for a given Mach number, we start with an initial guess 
for the equilibrium temperature. We then iterate by find- 
ing the equilibrium temperature achieved by matching the 
net gaseous phase molecular cooling ( |Goldsmith|200l"[) with 
net heating, including turbulent heating function (Pan & 



Table 1. Fiducial Scalings Used in this Paper 



Padoan 2009[ ), and dust grain-gas cooling (HoUenbach & 
McKee|1989[ ). During the calculation, we maintain the dust 
grains to be at 10 Kelvin, as dust grains can radiate away 
heat very efficiently at the GMC densities. An updated Mach 
number is computed at an isothermal sound speed reflective 
of the new mean temperature. We then iterate until we reach 
convergence, ending with a self-consistent choice for both the 
temperature and Mach number, which will in general both 
differ from their initially assumed values. The converged val- 
ues for both the mean temperature and the Mach number 
are then used for the Eulerian calculation. For the simu- 
lations presented here, we have employed a self-consistent 
mean temperature of 14.8 Kelvin at a 3D RMS Mach num- 
ber of 3.5. 



3.1.2 Scaling 

Our simulated models are cast in dimensionless units such 
that the problem domain size L = 1 in each of the three 
spatial dimensions, the isothermal sound speed Ciso = 1, 
and the total mass M = 1. In this formulation, only one 
dimensionless parameter - the 3D RMS Mach number A4. - 
fully specifies the hydrodynamics of the problem for a given 
turbulent realization. 

It is convenient for reference to scale the dimension- 
less problem to physical units using fiducial values. In order 
to facilitate this conversion, we utilize the constraint that 
GMCs are observed to be in approximate virial equilibrium 
between turbulent kinetic energy and gravitational binding 
energy. We define the virial parameter Oviriai of a spheri- 
cal cloud of mass M and radius R to be the ratio of its 
turbulent kinetic energy to gravitational binding e nergy : 
a^i.i^i = 5cr^i?V(2GM) ([Bertoldi fc McKee|[l992|, where 



(T = M Ciso/3 is the 1-D velocity dispersion. We general- 
ize the definition to a uniform periodic domain by setting 
L = 2R, 



3-D RMS Mach number M 


5 


Isothermal sound speed Ciso 


2.3 X 10"* cm/s 


3-D RMS velocity vq 


8.07 X 10* cm/s 


Domain size L 


0.405 pc 


Domain mass M 


102 AfQ 


Dynamical time L/vq 


0.43 Myr 



Once we require that the cloud be in approximate virial 
balance with Qviriai = 1, this dimensionless model can be 
rescaled to values of direct astrophysical relevance. We can 
perhaps most easily see this by fixing fiducial values for the 
isothermal sound speed Ciso and the mean cloud density po- 
The domain size L then becomes 



,-1^, (12) 

DGpoQvirial 

where po is the mean cloud density, and the total mass M = 
poL^. The 3-D RMS velocity Wo on the domain scale is then 
simply M Ciso- 

In this paper, we choose to scale to a mean density 
typical of GMC clumps, with number density 4 x 10'' cm~'^ 
and mass density po = 1.28 x 10^^^ gm cm~^. This density 
corresponds to mean column density 'Eq — poL — 0.16 gm 
cm~^. We summarize the fiducial scalings used in table [I] 



3.2 Driven Turbulence 

Because turbulence is an inherently dissipative phenomenon, 
we drive the simulation with momentum and energy forcing 
terms in order to achieve a statistically-steady state. The 
stochastic driving method used in our simulation is the same 



'^virial 



bATci^ 
6 GM 



(11) 



as originally proposed by Eswaran & Pope ( 1988 1. The tur- 
bulent velocity fluctuations are described by Fourier trans- 
form from the spatial domain. For each stirred mode of the 
velocity fleld, an acceleration is applied at each time step. 
The fleld consists of three complex phases, with each acceler- 
ation mode evolved by a Ornstein-Uhlenbeck (OU) random 
process which is analogous to Brownian motion in a viscous 
medium. An OU process is defined to be a time-correlated 
time sequence, and has the properties of having a zero mean, 
and a constant root-mean-square in time. Each step in the 
sequence begins with the previous value, adds a Gaussian 
random vari able wit h a given variance, weighted by a driv- 
ing factor y/l — f'^ where / = exp (— Af/rcoi ) then decays 
the previous value by a factor /. Here, the OU process repre- 
sents a velocity forcing, and its variance is chosen to be the 
square root of the specific energy input divided by the de- 
cay time Tcor. In the time limit that the time step At — >■ 0, 
the algorithm represents a forcing term which is a linear 
weighted summation of the old state with the new Gaussian 
random variable. 

At each timestep, the velocities are then converted to 
their physical space by a direct Fourier transform, adding 
the trigonometric terms explicitly. Since the stirring is done 
in the low-order modes, most driving involves a fairly small 
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number of modes. Therefore this decomposition is more ef- 
ficient than a fast Fourier transform. The range of modes 
over which turbulence is driven includes wavelengths 1 /4 to 
1 times the box size in the simulations presented here. In 
the turbulence literature, in which the box length is set to 
2n, this range of scales corresponds to wavenumbers k — 1 
to 4. We note that, contrary to previous claims, the acceler- 
ation field in the FLASH default turbulence driving module 
is non-solenoidal, with about a 10% admixture of power in 
compressional modes. 



3.3 Lagrangian Tracer Particles 

A simulation built upon the Eulerian framework of hydro- 
dynamics discretizes the volume of the fiuid into individual 
spatial cells with some finite width Aa;. The Eulerian frame 
necessarily artificially mixes fluid parcels with different his- 
tories with a resolution-dependent mixing process when they 
enter into the same cell. In contrast, a Lagrangian compu- 
tational framework differs from the Eulerian framework in 
a basic yet essential regard : individual fluid elements are 
followed as a function of time t, without any artificial mix- 
ing. The Lagrangian framework is therefore ideally suited 
for the purpose of following the chemical evolution of the 
fiuid without the artificial mixing imposed by the Eulerian 
mesh. Any mixing in the Lagrangian framework must be 
introduced explicitly, which greatly facilitates a controlled 
treatment of this complex process. 

Passive Lagrangian tracers are characterized by their 
kinematic data alone : their positions Xiit), and their veloc- 
ities Vi{t), where i is an index which runs over the parti- 
cle number, fully specify their evolution. In addition to this 
kinematic data, the particles also retain hydrodynamic data, 
including the density pit) captured from interpolation from 
the Eulerian mesh, for use in chemical post-processing. Sig- 
nificantly, Lagrangian tracers can be thought of as massless, 
passive particles. Consequently, while tracers are advected 
by the fiow, they do not directly infiuence the velocity or the 
pressure of the fluid, and therefore do not back-react upon 
the fiow. 

The particle positions Xiit) are governed by the time 
integration of their velocities Vi{t) obtained by interpolation 
from the mesh. The numerical implementation of the tracers 
therefore hinges directly upon both the interpolation scheme 
and the time-integration scheme, which we now detail. 



3.3.1 Particle Interpolation from Eulerian Mesh 

The software suite FLASH provides various interpolation 
methods for Lagrangian tracer particles to obtain hydrody- 



ing the profile of the scalar field within a mesh cell. For 
the Piecewise Parabolic Method (PPM) method, the recon- 



namic data from the Eulerian mesh Dubey et al. (20091; 



Fryxell et al. (20001. Volume-averaged Eulerian physical 



quantities are discretized onto the mesh in a cell-centered 
fashion. For example, a continuous scalar field f{x) in one 
dimension is discretized in a given cell i as: 



1 

Ax 



i+l/2 



f{xl)dxl. 



(13) 



where i is the cell index, and Aa; is the cell size. The cell 
extends over the interval from a;,_i/2 to Xi^-112. An approx- 
imation to the function f{x) is obtained after reconstruct- 



structed function is second-order accurate spatially ( Colella 
fc Woodward_1984) . In this work, we employ a quadratic in- 
terpolant scheme to achieve the same order of accuracy for 
Lagrangian tracer particles. 



3.3.2 Particle Time- Integration Scheme 

The time integration for passive particles in our simulation 
is implemented by Heun's method, a two-stage Runge-Kutta 
explicit integration scheme. For simplicity, again we present 
the algorithm in one spatial dimension. We represent the 
fiuid velocity at position x and time t by v{x,t). Because 
the Eulerian mesh retains only discretized information of the 
fiuid velocity, v{x,t) is obtained by interpolation from the 
mesh to the location of the particles as just discussed. We 
therefore implement Heun's method by integrating particle 
positions forward in time from time step n at time to 
time step n -I- 1 at time f""'"^ via 

^x- + ^[v-+vr+'] (14) 

Superscripts denote time steps; subscripts, as noted previ- 
ously, denote particle index. Here At" is the increment of 
time from the n step to the n + 1 step, At" = t""*"^ — f". 
The intermediate velocity v*'"^^ is given by 



(15) 



(16) 



*,n + l 

and the intermediate position x*'""^^ by 

* ,n+l n 1 A ,n n 

X- — Xi + At Vi . 

The final particle locations x""''^ and velocities w"^^ are 
stored at time t""*"^, along with hydrodynamic data inter- 
polated to the particle locations, for use in post-processing. 



3.4 Gas-Phase Chemistry 

The chemical reaction networks are typically stiff due to the 
great span of reaction rates involved. Furthermore, astro- 
chemical rates are often highly uncertain; if they are mea- 
sured at all, they are typically measured at room tempera- 
ture and then extrapolated down to the lower temperatures 
typical of GMCs. The physical constraints are combined 
with the need to flexibly generate networks of hundreds 
to thousands of reactions required to model the gas-phase 
chemistry of GMCs. To address these challenges in the gas 
phase chemistry, we have built a new astrochemistry code 
framework, ASTROCHEM, upon the chemical kinetics soft- 



ware suite called the Kinetic Preprocessor (KPP) (Damian 
et al.|2002||Sandu et al.|2003[|Daescu et al.|2003t . 

KPP was designed as a general chemical kinetic tool 
to facilitate the numerical solution of chemical reaction 
networks. KPP automatically generates FORTRAN code 
that computes the time-evolution of stiff networks of chem- 
ical networks, starting with a high-level specification of the 
chemical reaction and rates, and allows the user to select 
from a variety of stiff integration schemes, including the Liv- 
ermore Stiff Ordinary Differential Equation Solver (LSODE) 
package ( |Gear|1982| and Rosenbrock methods). Significantly, 
because it symbolically computes the Jacobian and Hessian 
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matrices of the networks needed to interface with numer- 
ical integration schemes, it is capable of both generating 
new network codes, and treating sensitivities automatically. 
Furthermore, KPP exploits sparsity in the Jacobian and 
Hessian matrices to maximize computational efficiency. The 
ASTROCHEM framework interfaces to existing astrochem- 



ical databases and codes, including UMIST (Woodall et al, 
20071) and Nahoon (IWakelam & Herbst||2008') . The result- 



ing framework features a high degree of flexibility in cre- 
ating astrochemical reaction networks and analyzing their 
corresponding sensitivities. In the present paper, we focus 
primarily on carbon, nitrogen, and oxygen ion-neutral and 
neutral-neutral reactions in the gaseous phase. 

The chemical evolution of any species in our chemical 
network frame can be written in mathematical form as: 



1 k 1 
X X 



(17) 



m=l7^i 



where in the symbol x\ X denotes chemical concentration 
of a species, superscript i denotes a particular species. The 
term X/j^i^^i ^rS^sXn represents formation of species x' from 
species x^ through one body dissociation reactions, with 
Ar being reaction rate coefficient for one body dissocia- 
tion reactions and ai is the associated stoichiometric co- 
efficient, and N is the total number of species. Similarly, 
X^k^i-^i X/iliT^i B'^'bsCsXnXn IS the creation term for forma- 
tion of species x' from species x** and xS with B''' being 
two body reaction rate coefficient, bg and c^ are associated 
stoichiometric coefficients. Reactions of these type may in- 
clude neutral-neutral, charge transfer, recombination reac- 
tions etc. We can account for the destruction reactions too, 
with D™ being reaction rate coefficient for two body de- 
struction reactions for species x\ E| is the reaction rate 
coefficient for single body reactions, f^ is the stoichiometric 
coefficient for single body destruction reactions, dl and e™ 
are the stoichiometric coefficient for two body destruction 
reactions. 



3.5 Gas-Phase Cooling 

We next describe the gas phase cooling function A{n, T) 
adopted in our models. We use nearly-identical cooling terms 



as those used by ( Smith & Rosen 2003 1. In our model, we in- 
clude gas-grain (dust) cooling pollenbach fc McKee|1989[ ), 
where we take the dust temperature to be Tdust ~ 10 K. 
We also include cooling due to coUisionally-excited vibra- 



tional and rotational modes of molecular hydrogen ( Lepp & 
Shull|1983| ), rotational modes of water ( [Neufeld fc Kaufman 



1993), and coUisionally-excited vibrational modes of wa- 
ter ( HoUenbach fc McKee|[T989 ) . We also incorporate cool- 
ing due to the dissociation of molecular hydrogen ( [Shapiro] 
fc Kangl 119871 iLepp fc Shulll 119831). Lastly, we include 



carbon monoxide cooling through several different chan- 
nels - including coUisionally-excited rotational modes , and 



collisionally-excited vibrational modes ( Neufeld fc Kaufman 



1993 Neufeld et al.|[T995 l. The dominant gas phase molec- 
ular cooling mechanism among those listed above is via CO 
(carbon monoxide) rotational modes excited by collisions 
with both atomic and molecular hydrogen, AcorotcoiH(M, T) 
( Neufeld fc Kaufman|p93| [Neufeld et al.||1995[ ). The cool- 
ing rate AcorotcoiH(/i, T) for temperatures below 100 Kelvin 
( Neufeld et al.|1995| , while the cooling rate AcorotcoiH(w, T) 



for temperatures above 100 Kelvin is defined by (Neufeld fc 



|Kaufman^^l993| ). Lastly, the gas phase cooling function also 
takes into account the negative of the hydrogen reformation 
heating iHoUenbach fc McKee[[1979 1. 



When calculating the CO cooling rate, we determine the 
logarithm of optical depth Top. The optical depth is defined 
by: 



loglO (^op) = logi 



nco 



(18) 



where nco is the number density of carbon monoxide 
molecules per cubic centimeter, and dv/dr is the velocity 
gradient. This is a local approximation which ignores the 
shielding effect of the cloud and ignores the internal struc- 
ture of the cloud. We use an average value of velocity gra- 
dient calculated over the entire mesh for the optical depth 
calculation ( Pan fc Padoan|2009 1 : 



dv 
dr 



dux dvy 

— - H 

dx dy 



dz 



(19) 



where the subscript mesh denotes that an averaging over the 
entire mesh has been done. 



3.6 Shock Modeling 

When post-processing the Lagrangian tracers, we utilize a 
one-dimensional model of radiative shock dynamics, includ- 
ing heating, cooling and gas phase chemistry ( [Hollenbach[ 
fc McKee|1979| [Smith fc Rosen|20"03| ). The time-dependent 
equations of hydrodynamics, which respectively represent 
conservation of mass, momentum, energy and molecular hy- 
drogen are : 



9p ^ djpv) ^ ^ 
dt dx ' 

d{pv) d{pv'^+p) 



dt 

djpe) 
dt 



+ 



+ 



dx 
d{pev) 
dx 



= 0, 



= -P 



dv 
dx 



and 



d{fn) , difnv) 



dt 



+ 



doL 



R{T,n,f)-D{T,n, /). 



(20) 
(21) 
(22) 

(23) 



Here n is the total number density of hydrogen nuclei, in- 
cluding both atomic and molecular forms of hydrogen, pe is 
the total internal energy density, / is the fractional molec- 
ular hydrogen number abundance relative to n, A is the 
specific cooling rate, R is the specific rate of formation of 
molecular hydrogen, and D is the specific rate of dissociation 
of molecular hydrogen. 

We can express the total number of particle per unit 
volume ntot in terms of the fractional helium abundance 
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with respect to hydrogen atomic nuclei /(He) and molecu- 
lar hydrogen abundance as ritot ~ n[l + /(He) -/]. Then, 
the ratio of specific heats 7 and the mean molecular mass 
jj., taking into account the atomic and molecular phases of 
hydrogen as well as the atomic helium, may be written as : 



7 = 



5.5-3/ 
3.3-/ ' 



ntotrrip 



1 + if (He) 
1 + f{He) - f 



(24) 



Here we have assumed that the gas kinetic temperatures are 
high enough to excite the rotational degrees of freedom of 
H2, but not high enough to excite the vibrational degrees of 
freedom. 

This time-dependent hydrodynamic shock model speci- 
fied by equations[20]-[23]are further simplified for solving the 
hydrodynamic flow in the post-shock region by assuming the 
shock is stationary. This assumption is valid for the purposes 
of post-processing because the flow timescale is typically 
much less than the CFL timestep on the mesh; equivalently, 
the post-shock cooling layer thickness is much less than the 
typical grid spacing. In steady state, the Rankine-Hugoniot 
jump conditions yield the hydrodynamic state immediately 
post-shock : 



(7 + l)X' 



Po 



VI 



(7- l)M2 + 2 



El 

PO 



l + il--)M 



To 



Pi 

Spo' 



(25) 



(26) 



Here S is the compression factor across the shock. Quanti- 
ties labelled with the subscript "0" represent the hydrody- 
namic state immediately pre-shock, and quantities labelled 
with the subscript "1" are immediately post-shock. Because 
we consider only non-destructive shocks in this paper, the 
molecular abundances remain continuous across the shock. 

Under the stationarity assumption, the equations of hy- 
drodynamics yield the conservation of mass 



pv = plVl, 

and momentum 

p + pv'^ =pi+ pivf, 



(27) 



(28) 



in a steady shock relate the density p, velocity v, and pres- 
sure p throughout the cooling layer to the immediate post- 
shock state. We may rewrite the total internal energy in 
favor of pressure, pe = p/(7 — 1), to express the energy 
equation [22] in the following form : 



dp/ (7 — 1) 7 9' 



+ 



1^ dx 



-A(T,n,/). 



(29) 



dx 7 

Applying the conservation of mass and momentum 



(eqns. 27 and 28 1 to equation 29 we obtain an ordinary 



differential equation for the rate of change of the ratio of 
specific heats 7 : 



97 
dx 



4.4 



df 



(3.3 -/)2 dx- 



(30) 



Similarly, applying conservation of mass and momentum to 
equation |23[ we obtain another ordinary differential equa- 
tion governing the composition fraction / of molecular hy- 
drogen : 



df _ pmp{R 
dx 



D) 



PlVl 



(31) 



This formulation yields two simultaneous ordinary dif- 



ferential equations (eqns. 30 and 31 1 governing the post- 



shock cooling layer, which we solve with the stiff ordinary 
differential solver LSODE. Our network has over 200 chemi- 
cal reactions and 66 species (seejAj. The electron concentra- 
tions, which represent the sole negatively-charged species in 
our network, are calculated by charge conservation assuming 
net neutrality. The cooling terms used are same as that used 
by Smith and Rosen (see[A|. We integrate forward through 
the shock until the post-shock temperature falls to within 
0.1% of the ambient temperature. We then match the post- 
shock integration condition back onto the Lagrangian tracer 
trajectory. 

At each timestep of the simulation, the Lagrangian 
tracer particle time history stores the vector position Xi{t) 
and the vector velocity Vi{t) of the ith particle in Cartesian 
coordinates, along with the Eulerian density field p{xi{t)) 
interpolated onto the location of the particle. Since we have 
assumed isothermal gasdynamics for the Eulerian calcula- 
tion, the Lagrangian tracer particles must cool back down 
to the ambient temperature in order to maintain consis- 
tency between the Eulerian and Lagrangian fields. Due to 

the highly-intermittent nature of supersonic isothermal tur- 
bulence, statistically we expect to find portions of the GMC 
where the density lies at the extremely rarefied low-end 
tail of the density probability distribution function (PDF) 



( Vazquez- Scmadcni " 1994; Padoan & Nordlund' 19971 IPassot 



& Vazquez-Semadeni^l998j ^ederrath et al. ^20081. At these 



low densities, the local cooling time scale of the gas can 
be greater than the global dynamical time. In gas parcels of 
molecular clouds of number density n ~ 1 cm~^, the phys- 
ical cooling time due to rotational line cooling of carbon 
monoxide colhsionally excited by atomic hydrogen is ~ 3 
million years, assuming /(CO) --^ 7 x 10~^. We estimate an 
upper-bound to the cooling time tcooi of the gas by the ratio 
of the excess internal energy of the gas above the ambient 
internal energy Ae to the specific cooling rate A(n, T), tak- 
ing the collisionally-excited rotational line cooling of carbon 
monoxide as the dominant cooling mechanism. : 



tc 



Ae [1.5n(He) + 2.5n(7/2)] x fcs(Ti - To) 



A{n,T) 



A 



COrotcolH 



(n,r) 



(32) 



When tcooi < shock dynamical time, we identify density 
jumps as isothermal J-shocks. When this condition is not 
satisfied, we integrate forward in time without requiring the 
temperatures to return to their initial isothermal value. 

When postprocessing the trajectory, we filter the den- 
sity using a boxcar average, to help easily distinguish true 
shocks, and eliminate small post-shock oscillations. The 
derivative of filtered density is calculated in a geometrically 
slope-limited way. A negative time derivative of density rep- 
resents compression. On finding a continuous region of com- 
pression, the maximum physical density inside the region 
and the density just before entering this region are found. 
We subject these densities to the Rankine-Hugoniot shock 
conditions to retrieve the strength of the shock, assuming 
this compression to the be effect of a shock. We classify a 
continuous region with a negative density time derivative 
and compression ratio greater than 2 as a physical shock. 
We have plotted figures [2] and [3] as an example of how the 
algorithm works. We take a Lagrangian trajectory, and plot 
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the number density variation as a function of time in fig- 
ure [2] Plotted in another panel in the same figure are the 
time instants where our algorithm detected a shock, and the 
vertical axis depicts the Mach number of the shock our algo- 
rithm detected. In figure[3] the pre-shock and the post-shock 
densities of the shocks have been plotted which were shown 
to be detected in figure [2] 

We implement a box-car average of the form p" = A x 
p""^ -|- (1 — A) X p", where p is the filtered density and p is 
the unfiltered density, and the superscript denotes the time 
step. Because a sharp discontinuity is diffused on the mesh 
after filtering, the filtering parameter A dictates the spatial 
resolution of shocks. We construct a filter which preserves 
well-separated distinct shocks and avoids artificially merging 
them. Our criterion for well-separated distinct shocks is that 
they must be separated by the inverse Nyquist wavenumber 
on the mesh, which we take to be 4 grid cells. The relation 
between the filtering parameter A and the number of cells 
A*' over which the discontinuity broadens is given by : 



A = (1 - fti 



(33) 



where f threshold is threshold level, close to unity, which iden- 
tifies the discontinuity. We set the smoothing parameter 
A by considering the action of the filter upon a Ifeaviside 
function. In all runs described here, we utilize a value of 
A = 0.3862, which returns a Heaviside function to a value of 
fthreshoid = 0.99 over four cells. 



3.7 Solenoidal Heating 

In the previous subsections, we have taken the compressional 
heating due to shocks into account. We now illustrate the 
steps taken to account for solenoidal heating due to dissipa- 
tion in vortex filaments and shear layers. 

The sum total of the spatially-and time-dependent 
solenoidal and compressional heating per unit mass e{x, t) 
on a mesh is represented as ( [Kritsuk et al.|2007[ ): 

e{x,t) = ^^^{\Vxu\^ + ^\V-u\'y (34) 

where Re is the Reynolds number, |V x u\^ represents the 
solenoidal kinetic energy, ||V ■ up represents the compres- 
sional kinetic energy. 

Because our simulations assume the inviscid Eulerian 
equations of hydrodynamics, we must infer an effective 
spatially-averaged Reynolds number Re which treats the nu- 
merical and artificial dissipation on the mesh as a real physi- 
cal viscosity. To do so, we average the dissipation over space, 
calculating the ratio of the spatially-averaged sum total of 
solenoidal and compressional kinetic energy by the average 
value of total dissipation e, 



Re = 



(iVxup + llV-ul^) 



(35) 



The mean dissipation e is ~ Vrms/iintograi. Using equation 
35 we find the averaged Reynolds number of a 512^ simu- 
lation of 3-D RMS Mach number 3.5 to be ~ 4200. 

An important point to consider is the relative strength 
of the solenoidal and compressional heating. We consider the 
ratio res of kinetic energy within compressible modes to the 



total kinetic energy, which is closely related to the ratio of 
compressible heating to total heating (eqn. [34|): 



res 



(|V-H^) 



(|V-u|2) + (|VXKj2)' 



(36) 



The ratio res varies between 0.28 (Kritsuk et al. 20071 to 
0.33 (Elmegreen & Scalo 2004a|b[ ) in simulations in the lit- 
erature. The mean value we obtain is 0.32, in close agree- 
ment with those previously reported. This ratio is subject to 
small-scale fiuctuations, and hence to resolution effects. An- 
other way to look at this is the ratio between compressive en- 



ergy and total energy ( Federrath et al.|20I0) , with compres- 
sive energy = 1(1^' w]"^), solenoidal energy = ^|V x m^^, 
and total energy being equal to sum of solenoidal and com- 
pressive energy. This ratio for our case is equal to 0.384. 

Significantly, the ratio res — 0.32 indicates that the ma- 
jority of the turbulent heating within the domain is due to 
the solenoidal dissipation. This result motivates us to treat 
the solenoidal dissipation on resolved scales only, using no 
subgrid modeling. While we expect the effects of solenoidal 
heating in our models will not be as pronounced as those 



works which perform subgrid modelling of vortices (Pan & 



Padoan 2009 Godard et al. 20091, the total heating rate 



is significantly higher than models treating shock heating 
alone. 

We now present the steps for calculating the evolution 
of the specific internal energy of the gas change outside of 
shocks. We adopt a Lagrangian equation for internal energy 
evolution outside of shocks, accounting for solenoidal dissi- 
pation : 



de 
di 



esix,t) + - (Tcrin) - A{n,T)) 
P 



(37) 



Here de / dt is the rate of change of specific internal energy, 
Es {x, t) is the solenoidal dissipation rate per unit mass at 
the particle's position x{t). Vciiji) is the volumetric heating 
rate due to cosmic ray heating, and A(n, T) is the volumetric 
gaseous cooling rate - all of which are converted from func- 
tions of number density n and temperature T to Lagrangian 
rates per unit mass by dividing by the mass density p. This 
internal energy equation is solved separately for each parti- 
cle using the LSODE stiff implicit solver outside of shocks, 
along its entire trajectory. 

The volumetric cosmic ray heating in the Eulerian 
frame due to the kinetic energy released in the ionization 
of helium and hydrogen molecules is given by: 



= 7?nC(l -I- ./(HE)) erg cm"^ s" 



(38) 



Here, r) is the kinetic energy released for every ionization 



electron released, which we take to be 20 eV (Goldsmith 
[2001), n is the number density, C, is the ionization rate, 
/(HE) is the relative abundance of helium nuclei compared 
to hydrogen nuclei. 

Using typical values, we scale the cosmic ray heating in 
the Lagrangian frame : 



c 



X 10- 



1.14 



erg g 



(39) 
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The symbols in equation [39] are same as defined earlier. 

Similarly, the Lagrangian heating rate due to hydrogen 
molecule reformation Tu^ reform /p is adapted from HoUen- 
bach and McKee (iHoUenbach & McKee||1979 1. 



For comparison, we estimate the solenoidal dissipation 
from the total turbulent dissipation rate e{x, t). The fiducial 
scaling of the total dissipation e used in our work is derived 
in equation |B5| for a turbulent RMS Mach 3.5 GMC clump, 
the mean rate is 2.02 x 10~*erg g~^s~^. Assuming that the 
solenoidal mode of dissipation is two thirds of the total dis- 
sipation, the fiducial value of the average physical solenoidal 
dissipation Is is then: 



= 1.34 X 10"*'erg g~^s~\ 



(40) 



On comparing equations |40] and |39| we find that the 
dissipation due to the solenoidal on the resolved scales with 
no subgrid modelling is of the same order of magnitude as 



the cosmic ray heating on average (Pan & Padoan 2009), 



though unlike cosmic ray heating, the turbulent dissipation 
is highly intermittent in space and time. Within the most 
intense structures of dissipation, the solenoidal dissipation 
rate es(x,t) will exceed the cosmic ray heating rate Fcr/p- 



4 VERIFICATION 

Our shock detection and chemical evolution algorithms de- 
pend on the ability of the Lagrangian tracer particles to 
accurately capture the hydrodynamical evolution. Conse- 
quently, the verification of the code modules plays a cru- 
cial role in understanding both the successes and the limi- 
tations of underlying numerical methods. The FLASH code 
is the product of nearly a decade of intensive software de- 
velopment, including careful attention to code verification 
(Rosner et al. 20001. FLASH is tested nightly against a 



suite of tests to verify the correctness of core physics mod- 
ules, including the PPM hydrodynamics module and the La- 
grangian tracer module. Standard nightly tests include the 
Sod shock tube problem [Sod^ (1978), the Colella- Woo dward 
colliding blast wave problem (Woodward & Colella 19841, 
and the Taylor-Sedov blast wave problem ( |Sedov||1946| ). In 
addition to these standard tests, we also constructed and 
performed an inclined strong adiabatic shock test to verify 
the accuracy of Lagrangian tracer particle advection in the 
presence of strong shocks inclined with respect to the mesh. 

In addition, the same Eulerian PPM hydrodynamics 
and Lagrangian tracer modules used here were employed 
in a large-scale computational study of weakly-compressible 
turbulence on a 1856^ Eulerian mesh, with 256^ Lagrangian 



tracer particles (Benzi et al. 2008 20101. A stringent 



comparison study of our computational methodology for 
turbulence was recently completed by five computational 
and three experimental groups studying fundamental La- 
grangian structure functions of turbulent flows ( |Arneodo| 
et al.|2008 |. The excellent agreement obtained in that cross- 
comparison between both multiple computational codes and 
experimental results represents both a stringent verification 
and validation of both the FLASH Eulerian PPM hydro- 
dynamics solver and the Lagrangian tracer module in the 
subsonic, fully turbulent regime. In a new verification test, 
detailed below, we consider a rigorous Lagrangian test for 



fully-developed, isothermal supersonic turbulence by requir- 
ing that in a statistically-steady state, the one-point proba- 
bility distribution functions of the Eulerian and Lagrangian 
density values must agree with one another. 

Lastly, we also conducted extensive tests verifying our 
gas phase chemical network solvers and cooling terms. We 
compared against previous results for a static background 



(Iglesias 19771 and single shock chemistry (Iglesias & Silk 



^1977fllsing ASTROCHEM. 



4.1 Density PDF Verification Test for Supersonic 
Turbulence 

Moving beyond verification tests in simple geometries, we 
next consider fully-developed supersonic turbulence. In or- 
der to verify our Lagrangian tracer methodology in the 
regime of fully-developed turbulence, we utilize a key iden- 
tity; namely, that the one-point Eulerian and Lagrangian 
PDFs in homogeneous, isotropic turbulence must be iden- 
tical ( Pope|[2b00 1 . Specifically, while both the Eulerian and 
Lagrangian PDFs derived from the simulation will in general 
differ from reality, they must both agree with one another 
by general consistency requirements. In addition, it is im- 
portant to note that this equivalence between the one-point 
PDFs holds rigorously, and provides the basis for a rare ex- 
act verification test in fully-developed turbulence. 

There are in fact two levels of verification involved in 
this test. In the first, exact level, we confirm that the La- 
grangian PDF derived from the simulation is statistically 
identical to the Eulerian PDF, also derived from the same 
simulation. Because the Lagrangian tracers follow the mass 
distribution, this test is a stringent verification of the La- 
grangian evolution. In particular, an arbitrary random dis- 
tribution of the tracers throughout the spatial domain will 
not in general succeed in reproducing the Eulerian distri- 
bution. In the second level of the test, we verify the Eule- 
rian PDF against the body of literature of numerical simula- 
tions, which have demonstrated that supersonic isothermal 
turbulence yields an Eulerian mass density probability dis- 
tribution function which follows a log-normal distribution 



(|Vazquez-Semadeni|1994 Padoan fc Nordlund|1997 Passot 



fc Vaz quez-Se madeiiil 19981 INordlund fc Pad oan'199 9| |Kr~ 
suk et al. 2007| [Lemaster fc Stone||2008| |Fcdcrrat h et al. 



2008 



2010) 



The Eulerian PDF is calculated by equally binning the 
density between the maximum and minimum values. For the 
Lagrangian PDF, we weight by mass by weighing the num- 
ber of particles falling into each density bin by 1/p, where p 
is the density of the Lagrangian tracer particle being added 



to the bin ( Kritsuk et al.^,2007p . We list the important run- 
time parameters of the simulation in table [2] All runs utilize 
a ratio of specific heats 7=1-1- 10~® and an Ornstein- 
Uhlenbeck turbulent driving correlation time tcorr = 0.5 
in dimensionless units such that the global sound-crossing 
time is unity. We allow the system to evolve over four dy- 
namical times to reach a steady state. We confirmed that 
the system has in fact reached a steady-state by verifying 
that the RMS velocity of the system approaches a near- 
constant value. Time-averages are then performed over the 
last two dynamical times. As a demonstration of our cal- 
culations where we convince ourselves that Lagrangian and 
Eulerian PDF come from the same underlying distribution, 
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Table 2. Results for Eulerian/Lagrangian PDF Verification. 



Table 3. List of chemical species evolved in the model 



Simulation Parameter 


Run 1 


Run 2 


Run 3 


7V3 . , 

grta 


128^ 


2563 


5123 


, . , 

partzcle 


643 


643 


643 


d-statistic 


0.0435 


0.0497 


0.0435 


KS probability 


0.9976 


0.9866 


0.9976 


b 


0.4988 


0.4551 


0.4351 



we have plotted the logarithm of Eulerian and Lagrangian 
PDF versus logarithm of density for a 256^ and 5123 gjjjj^_ 
lation in figure [4] 

We perform a best-fit to the time-averaged density 
probability distribution in log-normal form : 



p{ln p) d(\np) 



y(27ra2 



exp 



1 / In p — In p 
2 



d(lnp)(41) 



Here the mean of natural logarithm of density is given by 
h^=-^, (42) 



where p is the volume density ( Kritsuk et al. 2007 Federrath 
et al.|2008 20101. The relation between standard deviation 
o" and Mach number M for the parameter b is 



o-^ = In (1 + b^M^) 



(43) 



For our highest resolution simulation, we can obtain the 
parameters of the fit in two independent ways which are 
formally identical for an exactly log-normal density distri- 
bution. However, the intermittency inherent in turbulence 
causes the tails of the distribution to differ from an exact 
log- normal ( [Federrath et al.|2008||2010[ ). In the first param- 
eter determination, we simply determine the mean density 
directly to find In p = —0.9227. From equation 42 and equa- 
tion [43j we then obtain a — 1.3584 and b = 0.6596, in good 
agreement with previous authors. Alternatively, if we fit the 
density distribution to equation|4l]to obtain a = 1.3896, we 
find b — 0.6938 from equation |43[ We would like to point out 
that value of parameter & ~ 1 for purely compressive forc- 
ing and 6 ~ 1/3 for a purely solenoidal forcing (Federrath 



|et al.|[2008[ |2010"1 ). This apparent conflict in the determined 
values of b points to intermittency inherent in turbulence 
( Federrath et al.|2008[|2010[ ). 

Both the Eulerian and Lagrangian PDF distributions 
we obtain show significant departures from a log-normal 
fit at low densities. While early low-resolution studies pro- 
duced a log-normal fit to the density PDFs centered around 
the mean, more recent simulations at higher resolution have 



tended to be skewed to lower density values (Padoan & 
|Nordlund|[T997| [Kritsuk et "al]|2007| ). The low-density tail 
is more subject to turbulent intermittency, and so we ex- 
pect that the actual PDF should skew to low density val- 
ues. Kritsuk et al 2007 have conducted perhaps the highest- 
resolution study of isothermal supersonic turbulence statis- 
tics at = 6 on a 2048^ mesh, averaged over many snap- 
shots, accumulating the PDF for lO^'^ data points. The 
skewness seen in Kritsuk et al's results has since been con- 
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firmed by many other grid and SPH codes (Kitsionas et al. 
'20091 . Moreover, the trend is for higher resolutions to be 
even less skewed towards low density in the mean, because 
of the greater sampling of the intermittent tails of the turbu- 
lent distributions, which is the case for both grid and SPH 
based codes ( Price fc Federrat"li|[2010 1. 

We use the Kolmogorov-Smirnov(KS) test to check for 
the probability that both the Lagrangian and Eulerian prob- 
ability distribution are drawn from the same underlying den- 
sity distribution. The KS test is based upon the d-statistic, 
which is defined as the maximum difference in the cumu- 
lative distribution functions of the probability distributions 
being compared. As a consequence, the KS test is insensi- 
tive to any binning of the underlying PDFs, as well as any 
presumed functional fit. At our highest resolutions of 5123, 
we find a d-statistic of 0.0435, which implies a KS probabil- 
ity that the Lagrangian and Eulerian density PDF originate 
from the same underlying distribution with probability of 
0.9976. This result is a strong, rigorous verification of both 
the Eulerian and the Lagrangian hydrodynamics in a super- 
sonic isothermal turbulent fiow. 



5 RESULTS 

5.1 Time-Evolution of a Single Lagrangian Parcel 

We present the detailed chemical evolution along one tracer 
particle trajectory, under conditions typical of GMC clumps. 
The initial state is in chemical equilibrium at the back- 
ground state. We begin the simulation in chemical equilib- 
rium on the uniform background state of the cloud. Changes 
in chemical concentrations come about due to the density 
fiuctuations in the background cloud, as well as the tem- 
perature enhancements in s which a gas parcel encounters. 
The tracer particle is taken from a 5123 j-ggQjution Eulerian 
simulation, with an adiabatic coefficient of 7 = 1-1- lO"'', 
in a supersonic turbulent medium, with a 3D RMS Mach 
number of 3.5. We post-process Lagrangian trajectories as 
described in section |3.6[ in calculating the chemical evolu- 
tion for our chemical network. The complete list of chemical 
species evolved in our simulation are listed in the table [3] 
The time evolution of the relative concentration of some 
species along a single Lagrangian tracer are shown in fig- 
ure [5] The initial chemical concentrations, relative to the 
total number of hydrogen nuclei, used in deriving an initial 
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chemical equilibrium at a fixed background number density 
of 2 X IC* cm~^ and ambient interstellar temperature of 10 
Kelvin are listed in table lAlT] 

Inspection of figure [s] reveals that some species (e.g., 
CH2, HCO"'') experience rapid changes in relative chemical 
concentrations, in some cases of up to three orders of mag- 
nitude compared to their relative concentration in a static 
medium. Other species (e.g., H2O and NH3) evolve more 
smoothly, over a dynamic timescale. Still others (e.g., CO) 
remain relatively constant. 

To better understand the chemical evolution of these 
species, we identify two characteristic time scales, which 
jointly determine the chemical and thermodynamic evolu- 
tion in the post-shock cooling layer. The first time scale is 
the reaction time, which is defined for a given species X to 
be the inverse of the the physical rate of change of concen- 
tration of the species, normalized to its own number density 
: tchem(^) = nx/ (dnx/dt). Here nx, or [X] is the number 
density of a given species X. The second time scale is the 
cooling time defined as : fcooi = T/{dT/dt), the time in 
which gas phase temperature reaches the background tem- 
perature of the molecular cloud. 

When computing the cooling time, we have taken the 
dominant cooling rate to be that due to CO rotational 
lines excited by collisions with atomic and molecular hy- 
drogen. Assuming molecular hydrogen and helium to be 
the dominant species, the internal energy required for a 
parcel of gas in the immediate post-shock temperature 
Tp to drop down to the ambient temperature Tambiont is 
(1.5[He] + 2.5[H2]) x k^iT^ - Tanibicnt)), while the cooling 
rate is given by A(n, T) (see section 3.5 1. The cooling time 
in the post-shock flow is then given by 



cooling — 



(1.5[He] -I- 2.5[H2]) x k^ij^ - T^mhicnt)) 



(44) 



A(n,T) 

This cooling time is the same for all the chemical species. 

The chemical time tchom and the cooling time tcooUng 
can be used to define a local Damkohler number Da(X) in 
the post-shock fiow for a given species X as 

^cooling 



Da(X) 



(45) 



^chem (X) 

Significantly, the local post-shock Damkohler number 
of a given species specifies whether the chemistry of that 
species is either fast or slow, relative to the cooling time in 
the post-shock flow. Speciflcally, for Da >> 1, the chemical 
timescale tchcm (X) is rapid in comparison to the post-shock 
cooling time tcooiing. The chemical evolution of such large 
Da species is highly sensitive to the presence of shocks. In 
contrast, for Da << 1, the chemical timescale is slow in com- 
parison to the cooling time. In effect, these low Da species do 
not "see" the shocks, and consequently, the chemical evolu- 
tion of these species are insensitive to the presence of shocks. 

We now directly compare three representative species - 
methylene radical (CH2), water (H2O), and carbon monox- 
ide (CO). We calculate the cooling time and chemical time 
associated with each shock for the sampled trajectory. In 
flgure[6] we show results for methylene radical CII2. Figure 
[6] shows that the Damkohler number of methylene radical 
Da(CH2) lies in the range 1.24 < log(Da(CH2)) < 4.01. The 
initial chemical concentration of methylene radical relative 



Table 4. Number Abundances Relative to Hydrogen Nuclei 


Element 


Relative Number Abundance 


Atomic hydrogen H ... 


0.4 


Molecular hydrogen H2 


0.3 


Helium He 


0.14 


Atomic oxygen O 


1.17 X 10-" 


Atomic carbon C 


1.47 X lO-'^ 


Atomic sodium Na .... 


1.02 X 10-'' 


Atomic nitrogen N .... 


2.14 X 10-5 


Carbon monoxide CO . 


5.87 X 10-5 



to the total number of hydrogen nuclei is 1.97x 10~^^. As the 
trajectory evolves, the maximum relative chemical concen- 
tration of CII2 is 2.28 X 10"^^, while the minimum relative 
concentration of methylene radical is 1.64 x 10~^*. We infer 
that the relative concentration of CH2 increases by a max- 
imum factor of ~ 11.5 when compared to the initial static 
equilibrium state, while the relative concentration of CH2 
decreases by up to a factor of ~ 1.2 x lO'^x when com- 
pared to the initial static equilibrium. The time-averaged 
Damkohler number of methylene radical within the imme- 
diate post-shock cooling layers is logio(Da(CIl2)) = 2.76. 

In contrast to the large Damkohler number of CII2 in 
flgure |6j we show the smaller Damkohler number of water 
Da(H20) in figure |6] We determine that Da(Il20) lies in 
the range —0.636 < log(Da(H20)) < 0.556 on a logarithmic 
scale. We also see that the Damkohler number associated 
with any shock for methylene radical is greater than that of 
water. This means that for any shock, the equilibration time 
for water is much larger than that of methylene radical. The 
initial chemical concentration of water relative to the total 
number of hydrogen nuclei is 3.94 x lO"'^. As the calculation 
evolves, the relative chemical concentration of H2O attains 
a maximum of 4.80 x lO"'^. The time-averaged Damkohler 
number of water calculated within the immediate post-shock 
layers is logio(Da(H20)) = -0.02. 

Lastly, some species like CO (figure [5|, have nearly- 
constant relative chemical concentrations as the simula- 
tion evolves. As an illustrative example of this class of 
species, we consider CO here. From figure [6] we de- 
termine that the Damkohler number of carbon monox- 
ide Da(CO) on a logarithmic scale lies in the range 
-3.02 < log(Da(CO)) < -1.39. The relative concentration 
of carbon monoxide varies between 5.81 x 10"^ and 5.85 x 
10~^, which translates to a fractional change of ~ 6.8 x 10""^ 
in the relative concentration of carbon monoxide as com- 
pared to its static equilibrium state. The time-averaged 
Damkohler number of carbon monoxide within the immedi- 
ate post-shock cooling layers is logio(Da(CO)) = —1.99. We 
note that our results for the behavior of carbon monoxide 
differ from that of ( [Glover et al.|2010| ), as we only consider 
gaseous phase reactions, excluding external surface chem- 
istry and FUV radiation. The specific conclusions for CO 
hold for this class of species. Specifically, for species whose 
chemical time scale is disproportionately larger than the 
cooling time, there is virtually no change in the relative con- 
centration of the species throughout its evolution 

Based on these results, we classify the evolution of 
chemical species in a supersonic turbulent flow according 



14 A. Kumar and R. Fisher 



Table 5. Chemical Regimes Divided by Damkohler number 



Damkohler number 


Physical Significance 


0.50 ^ logioDa 


Shock-enhanced, stochastic variation 


-1.50 logioDa s£ 0.50 


Slowly-varying species 


logioDa ^ -1.50 


Frozcn-in abundance 



to their Damkohler number. In much the same way as vari- 
ations in dimensionless numbers characterize the transition 
between physical regimes in other hydrodynamic flow prob- 
lems - the most notable of course being the transition to 
fully-developed turbulence characterized by the Reynolds 
numbers - we characterize the transition in the character of 
post-shock chemical evolution of a species by its Damkohler 
number. We list our empirical findings in table |5] where 
we identify whether the concentration of a species is either 
shock-enhanced, slowly-varying, or frozen-in, depending on 
its Damkohler number. The ranges are drawn from the re- 
sults in this section, and are intended to be rough quantita- 
tive guides to the evolutionary regimes. 



5.2 Time-Evolution of an Ensemble of Lagrangian 
Parcels 

To gather a fuller picture as to the range of variation present 
in the physical and chemical structure of the GMC clump 
analyzed in section [STlj in this section we present the aggre- 
gate statistics for an ensemble of 100 post-processed parti- 
cles. First, in order to analyze the thermodynamic structure 
of the clump, we construct a joint PDF of density and tem- 
perature f{p, T), defined such that the integral over all tem- 
peratures of the joint PDF returns the PDF of the density 
field : 

poo 

fip)= / fiP,T)dT (46) 
Jo 

Conversely, integration over density returns the PDF of the 
temperature field f{T). 

In order to compute the Eulerian joint probability dis- 
tribution function of density and temperature shown in the 
figure [7] we bin the trajectories in a 500 x 500 matrix of log 
density and log temperature, weighing densities by a factor 
of 1/p. The entire joint PDF is accumulated by a random 
sample of 200 particles, averaged over all times. 

Figure [7] depicts the joint probability distribution func- 
tion f{p,T). Significantly, because we allow both time- 
dependent heating and cooling of the gas, the gas temper- 
ature is a function of the history of the fluid parcels, and 
therefore the joint PDF is not distributed along an equa- 
tion of state prescribed by a single adiabat , as it is in sim- 
pler models derived from time-independent thermodynamic 
equilibrium models of the gas I jSpaans fc Silk|2000| ). 
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Figure 1. A diagram illustrating the large dynamic range between the global dynamical time and the local thermodynamic and chemical 
time in post-shock cooling flows, for a representative species (CH2). The figure at upper left shows the evolution of the normalized 
abundance of CH2 over a global dynamical time, along a single fluid parcel. The vertical lines indicate the portion of the time domain 
blown-up at right, depicting the temperature and gaseous number density evolution post-shock. The flgure at lower left depicts the log 
normalized abundance evolution of CH2 over this same shock. 
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Figure 2. Shock location algorithm applied to an individual Lagrangian tracer particle. The top plot shows gaseous number density 
(in cm~^) versus time. Local Mach number versus time is depicted in the bottom plot. 
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Figure 3. Shock location algoritiim applied to an individual Lagrangian tracer particle. The top plot shows pre-shock gaseous number 
density (in cm~^) versus time. The bottom plot shows post-shock gaseous number density (in cm~^) versus time. 




Figure 4. Comparison of Lagrangian and Eulerian probability distribution function (PDF) for a 3D-RMS Mach number 3.5, supersonic, 
isothermal, turbulent flow. 
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Figure 5. Variation of normalized abundances of various species, relative to the total hydrogen nuclei density (= uh + nnj)) shown 
versus time for an individual Lagrangian tracer particle. Part 1 
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Figure 6. The top plot shows the variation in normalized H2O, CO and CH2 abundance versus time for a single Lagrangian tracer, in 
the presence of multiple shocks. The bottom plot shows the logarithm of local Damkohler number Da versus time. See text for definitions. 
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The joint probability distribution function is a power- 
ful tool which can be used to analyze a number of other 
questions, including the filling fraction of the warm molec- 
ular gas. By integrating over all densities, we obtain the 
temperature probability distribution function f{T). Using 
the temperature PDF, we find that 2.52% by volume of the 
molecular gas is at temperature of greater than 100 Kelvin. 
The temperature-integrated density PDF /(p) is approxi- 
mately log-normal. These results are in general accord with 
the established picture of shock-generated density fluctua- 
tions in supersonic turbulence ( Padoan fc Nordlund|[l997 1 . 
However, using our new methodology, we are able to more 
accurately quantify the temperature structure within thin 
shock and shear dissipation structures within the clump, a 
key step forward towards a more realistic chemical model of 
GMCs. 



6 DISCUSSION AND CONCLUSIONS 

Significantly, a key implication of this work is that the char- 
acter of the astrochemical evolution of a given species in 
a fully turbulent medium is robustly determined by its ef- 
fective Damkohler number, which can be simply estimated, 
without the need for the complex three-dimensional hydro- 
dynamic simulations, or even detailed shock calculations. 
This implication is of significance to observers. For instance, 
tracers like H2O, NH3, and Na^, which are less sensitive 
to the local turbulent fluctuations within a GMC, and are 
therefore more reliable tracers of the mean physical state of 
the cloud itself. Conversely, species like CH2, HCO'*', and 
OH are very sensitive to turbulent density fluctuations, and 
may serve as probes of the dynamics of turbulent intermit- 
tency. 

These results suggest that an adaptive reduced chem- 
istry model may be feasible in streamlining the compu- 
tational expense of a large chemistry network in three- 
dimensional hydrodynamical and magnetohydrodynamical 
models. In particular, by estimating both the local cooling 
time as well as the reaction timescale for a given species, one 
can determine the Damkohler number Da of that species. If 
Da >> 1, then that species can be locally approximated as 
being in chemical equilibrium. Similarly, species for which 
Da << 1 can be held fixed throughout the post-shock cool- 
ing layer. 

While in this paper we have treated only hydrodynamic 
J-shocks, we plan to extend this model to include magnetic 
turbulence, which will enable us to take into account the 
effects of partial ionization, including ambipolar diffusion. 
We note that while the normalized electron abundance in 
these models does not remain constant, it variation is lim- 
ited to about a factor of 2, from ~ 0.5 — 1.2 x lO"'^ (table 



All I. This suggests that a relatively straightforward exten- 



sion of this hydrodynamic methodology to non-ideal magne- 
tohydrodynamics by post-processing the chemical evolution 
along Lagrangian tracers may be possible. 

The non-equilibrium nature of the chemistry of some 
species may have particular significance for observational 
and theoretical studies of molecular abundances within star- 
forming GMC cores, and ultimately protostellar disks and 



planets. In particular, even in the absence of active star for- 
mation and outflows, the background turbulence within the 
GMC leads to significant departures of the abundances of 
some species from the chemical equilibrium at the mean 
cloud density. This variance in initial concentrations may 
have a direct influence on the zeroing of chemical clock 
timescales, which are typically taken about a chemical equi- 
librium state (Bcrgi n et al.|1997[ |Bcrgin fc Tafalla 2007') . In 
the context of a fully turbulent origin of GMC cores, a back- 
ground equilibrium state no longer exists, and the initial core 
abundances are themselves stochastically established by the 
shock dynamics. This effect may be particularly pronounced, 
if in fact GMC cores are formed via shock compression, as 
some leading models suggest ( Padoan fc Nordlund|1997 l. 
This stochastic variance of chemical abundances in the 



presence of shocks is in fact hinted at in earlier work by ( Ed- 
win A. Bcrginl & Ncufcld 1998) who studied the evolution 



of O2 and H2O under the influence of randomly sequenced 
C-shocks, as a model of stellar outflows, using Monte Carlo 
simulations. In their gas-phase models with signiflcantly 
stronger shock strengths than used here (up to 70 km s~^), 
they found a variation of a factor of ~ 3 in H2O abundances, 
and a barely noticeable effect in O2 abundances. Their find- 
ings are, however, consistent with ours in that we find both 
H2O and O2 to be low Damkohler number species at typi- 
cal background conditions within J-shocks, and also do not 
observe large variations in their abundances. Our results in- 
dicate that other species may be more sensitive to shock 
dynamics, even at the level of the background turbulence in 
regions away from active stellar outfiows. 
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0959382 and AFOSRDURIP Grant FA9550-10-1-0354. This 
research was supported in part by the National Science 
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Louisiana Optical Network Initiative under grant number 
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Table Al. Ionization reactions 



no. 




Chemical reaction 


rate coefEcient (s ^) 


Reference 


1) 


H2 - 


- c.r.p — >■ H^ + e~ + c.r.p/ .... 


0.95 X 10-1^ 


1 


2) 


Ha - 


- c.r.p — >■ H + H+ + e~ + c.r.p/ 


0.05 X lO-i'^ 


1 


3) 


He H 


- c.r.p — > He+ + e~ + c.r.p/ .... 


1 X W-^'^ 


1 



Table A2. Radiative Association 



no. 


Chemical reaction 


rate coefiicicnt(cm^s ^) Reference 


4) 


C+ + H2 CH2 + + hiy 


0.65 X 10-1^ 2 



^ Reference 1 jlglesias fc Silk|l977| l &; Reference 2 jWoodall et al. 
20071 
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Table A3. Positive Ion Molecule reactions reactions part 1 



no. 




Chemical reaction 


rate coefRcient(10 ^1 


cm^s ^) 


Reference 


5) 


Ha + 


H+ -5> H+ + H 


2.08 




2 


6) 


Ha + 


CO+ -5> HCO+ + H 


0.75 




2 


7) 


Ha + 


Na+ -5> HN2+ + H 


2.00 




2 


8) 


Ha + 


He+ ^ H+ + H +He 


3.7 X 10-5 X exp(- 


-35/T) 


2 


9) 


Ha + 


0+ -> OH+ + H 


1.7 




2 


10) 


Ha + 


OH+ HaO+ + H 


1.01 




2 


11) 


Ha + 


HaO+ ^ H3O+ + H 


0.64 




2 


12) 


Ha + 


N+ -5> NH+ + H 


exp {-85/T) 




2 


13) 


Ha + 


NH+ NHa+ + H 


1.28 




2 


14) 


Ha + 


NHa+ -5> NH3+ + H 


0.27 




2 


15) 


Ha + 


CH+ -5> CHa+ + H 


1.20 




2 


16) 


Ha + 


CHa+ CH3+ + H 


1.60 




2 


17) 


Ha + 


CN+ HCN+ + H 


1.00 




2 


18) 


Ha + 


HCN+ HaCN+ + H 


0.9 




2 


19) 


Ha + 


CNO+ -)■ HNCO+ + H 


1.0 




1 


20) 


Ha + 


HNCO+ HaNCO+ + H . . . 


1.0 




1 


21) 


Hc+ - 


f CO ^ C+ + + He 


1.6 




2 


22) 


Hc+ - 


f Na ^ N+ + N + He 


0.96 




2 


23) 


He+ - 


f Na -5> Na+ + He 


0.64 




2 


24) 


He+ - 


f O2 ^> 0+ + + He 


1.0 




2 


25) 


Hc+ - 


f O2 ^ Oa+ + He 


3.3 X 10-2 




2 


26) 


He+ - 


f CN ^ C+ + N + He 


0.88 X (T/300)- 


-0.5 


2 


27) 


He+ - 


f NO N+ + + He 


1.4 




2 


28) 


He+ - 


f NCO NCO+ + He 


1.0 




2 


29) 


He+ - 


f NCO ^ CO+ + N + He ... 


3.0 




2 


30) 


He+ - 


f NCO ^ N+ + CO + He ... 


3.0 




2 


31) 


He+ - 


f HNCO -5- NCO+ + H + He 


1.0 




1 


32) 


Hc+ - 


f HNCO HNCO+ + He ... 


1.0 




1 


33) 


He+ - 


f C02 ^ CO3+ + He 


0.121 




2 


34) 


He+ - 


f C02 -5> CO+ + + He ... 


0.87 




2 


35) 


He+ - 


f C02 ^ Oa + C+ + He .... 


0.04 




2 


36) 


He+ - 


f C02 ^ CO + 0+ + He ... 


0.1 




2 


37) 


He+ - 


f C02 -> Oa+ + C + He .... 


0.011 




2 


38) 


C+ + Oa ^. CO+ + O 


0.38 




2 


39) 


C+ + 


OH H+ + CO 


2.00 




1 


40) 


c+ + 


NH ^ H+ + CO 


2.00 




1 


41) 


c+ + 


CH -5> C2+ + H 


0.38 X (T/300)- 


-0.5 


2 


42) 


c+ + 


HaO HCO+ + H 


0.9 X (r/300)- 


0.5 


2 


43) 


c+ + 


NHa ^> H+ + HON 


2.0 




1 


44) 


C+ + CHa ^ CaH+ + H 


2.0 




1 
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Table A4. Positive Ion Molecule reactions reactions part 2 


no. 


Chemical reaction 


rate coe£Eicient(10 ^ 


cm^s ^) 




40 ) 


O ^ + xlL^iN — r Kj2^^ ~r ti 


0.1 X [1 / oUU ) 


-0.5 


Z 


40 ) 




1 1 
i. i 




z 


47) 


r~'+ 1 1\TU XT (^1^7+ 1 XT 


1 '7c; \y ('^^/Q^^^- 
i. /O X [I j oUU ) 


-0.6 


z 


4o ) 


' + rl2'^LJ — ? rL02vJ ' + xl 


n 




-1 
± 


^i) J 


nv^w ' wxi — r n\y\j2 -r n . . . . 


i.uu X y± /oxJKj) 


-0.6 


9 




■Lj/~ir\+ 1 "Lj V TJ r\-\- 1 r~*r\ 

ti\^\J ' -j- rl2*J —T rl3*J + ^--O . . . 


z.O X (^J / oUU ) 


-0.5 


z 


c;i ^ 
51j 


n \_ Kj ' -\- L\ri — r i>rl2 ~r W .... 


U.D4 X I^J /oUU J 


-0.6 


z 


OZ ) 


rlovj -h IN rl2 — 7' -lNrl3 -|- OL.* • . • 


u.oy X [I /ouu J 


-0.6 


z 


oo ) 


xlUvJ ' -|- lNrl3 — > lNrl4 ' -j- K^Kj . . . 


f I'' /QnnA — 
z.z X [I / oUU J 


-0.5 


z 


04 j 


tTr^o+ 1 XT (~*r^ \ XT i~^r\-\- _l /^/^ 
xlOvJ -|- rl3L^U ' + KjKJ 


0.0 X i^i / oUU J 


-0.5 


z 


00 j 


XT + 1 i^XT+ 1 XT 


U.o4 




z 


OO ) 


XT + 1 ^ r~'XT+ 1 XT 

-rl3^ + U ^ L^H^ + H2 






-1 
± 


57j 


XT + 1 r V ~» V xT/^r~»+ 1 XT 


i. / 




z 


ooj 


XT + 1 M \ XTM + 1 XT 

113^ + iN2 rliN2 + II2 


1.8 




2 


Oi)) 


XT + 1 ^^XT V XT /^+ 1 XT 


1.0 X [I /oUUj 


-0.5 


2 


fifW 
OKJ) 


XT + 1 (^M XT/^'\T+ 1 XT 

rl3 ' + L/iN — r rH^iN ' + rl2 


z.U X (^J / oUU ) 


-0.5 


z 


61) 


XT + 1 XT V XT 1 XT 

H3^ + II2U — )■ ll3(J^ + rl2 


0.9U X (-/ /oOUJ 


-0.6 


2 


62) 


H3+ + CO2 ^ HCO2+ + H2 .... 


2.0 




2 


63) 


H3+ + NH3 NH4+ + H2 


4.39 X (T/300)- 


-0.6 


2 


64) 


H3+ + HCN ^ H2CN+ + H2 .... 


8.10 X (T/300)- 


-0.5 


2 


65) 


H3+ + HCO ^ H2CO+ + H2 .... 


1.7 X (T/300)- 


-0.5 


2 


66) 


H3+ + H2CO ^ H3CO+ + H2 . . 


6.3 X (T/300)- 


-0.5 


2 


67) 


H3O+ + C -J- HCO+ + H2 


1.0 X 10-2 




2 


68) 


H3O+ + NH3 -J. NH4+ + H2O . . 


2.2 X (T/300)- 


-0.5 


2 


69) 


H3O+ + HCO H2CO+ + H2O 


1.0 




1 


70) 


O2+ + N ^ NO+ + 


0.18 




2 


71) 


HN2+ + CO ^ HCO+ + N2 


0.88 




2 


72) 


CH3+ + ^ H2CO+ + H 


4 X 10-2 




2 


73) 


CH3+ + N ^ H2CN+ + H 


6.7 X 10-2 




2 


74) 


il + ('{)■> JiCO - 






2 
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Table A5. Charge- Transfer Reactions 



no. 


Chemical reaction 


rate coefficient(10 ® 


cm^s ^) 


Reference 


75) 


H+ + O ^ 0+ + H 


0.731 X (r/300)0-23 exp 


-(225.9/T) 


2 


76) 


H+ + O2 ^ O2+ + H 


2.00 




2 


77) 


H+ + NO ^ NO+ + H 


2.9 




2 


78) 


H+ + OH ^ OH+ + H 


2.10 X (T/300)- 


-0.5 


2 


79) 


H+ + H2O H2O+ + H 


6.90 X (T/300)- 


-0.5 


2 


80) 


H+ + NH3 ->■ NH3+ + H 


3.70 X (T/300)- 


-0.5 


2 


81) 


H+ + H2CO ^ H2CO+ + H 


2.96 X (T/300)- 


-0.5 


2 


82) 


H+ + CN -5> CN+ + H 


1.0 




1 


83) 


H2+ + N -5> N+ + H2 


1.0 




1 


84) 


C+ + NO ^ NO+ + C 


0.52 




2 


85) 


C+ + NH3 -5> NH3+ + C 


2.2 X (T/300)- 


0.5 


1 


86) 


C+ + CH3 ^ CH3+ + C 


1.00 




1 


87) 


O2+ + C ^- C+ + O2 


5.2 X 10-2 




2 


88) 


O2+ + Na Na+ + 02 


0.71 




2 


89) 


N2+ + Na Na+ + N2 


2.00 




2 


90) 


NO+ + Na ^ Na+ + NO 


7.7 X 10-2 




2 


91) 


H2O+ + Na ^ Na+ + H2O 


6.20 




2 


92) 


H3+ + Na ^- Na+ + H2 + H 


2.1 




2 


93) 


H3O+ + Na ^ Na+ + H2O + H . . . 


3.1 




2 


94) 


HCO+ + Na ^ Na+ + CO + H ... 


2.60 




2 


95) 


HN2+ + Na Na+ + N2 + H 


0.50 




1 


96) 


NH3+ + Na ^ Na+ + NH3 


3.2 




2 


97) 


H2CO+ + Na Na+ + H2CO .... 


2.60 




2 


98) 


H3CO+ + Na Na+ + H2CO + H 


2.60 




2 


99) 


CH3+ + Na Na+ + CH3 


3.40 




2 


100) 


H2CN+ + Na Na+ + HON + H . 


0.5 




1 


101) 


CN+ + Na Na+ + C + N 


0.50 




1 
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Table A6. Ion-Electron Recombination Process 





Chemical reaction 






102) 


H+ + e" - 


^ H + hu 


3.5 X 10-''''(T/300)-°-'^5 


2 


103) 


Na+ + ■ 


Na. + hu 


2.7 X 10-"(T/300)-°-69 


2 


104) 


O2 + + 


-^0 + 


1.95 X (T/300)""-™ 


2 


105) 


NO+ + e" 


-s>N-|-0 


4.30 X (T/300)-"-'^^ 


2 


106) 


H3 + + e" 


-s> H2 + H 


0.234 X (T/300)""-'''2 


2 


107) 


H3 + + e- 


^H+H+H 


0.436 X (T/300)-'^-52 


2 


108) 


HCO+ + e- 


- CO + H 


2.4 X (T/300)-°-SS 


2 


109) 


HN2 + + e 


- N2 + H 


3.6 X {r/300)-°-5i 


2 


110) 


H3O+ + e" 


- H2O + H 


1.08 X (r/300)-''-s° 


2 


111) 


H3O+ + e" 


- _5.0H-|-H-|-H 


2.58 X (T/300)^°-50 


2 


112) 


NH3+ + e- 


- -s> N2 + H 


1.55 X (T/300)~°-5° 


2 


113) 


NH4 + + e 


- ^ NH3 + H 


13.69 X (T/300)-°-5 


1 


114) 


NH4 + + e 


- NH2 + 2H 


3.19 X (T/300)"''-''^ 


2 


115) 


H2CN + + 


e- HCN + H 


10.0 X (T/300)"''-5° 


1 


116) 


H2CN + + 


e- ^ CN + H + H 


10.0 X (T/300)-°-50 


1 


117) 


CH2 + + e 


- CH + H 


1.60 X (T/300)-°-60 


2 


118) 


CHs + + e 


- CH2 + H 


0.775 X (T/300)-°-80 


2 


119) 


CH3 + + e 


- ^ CH + H + H 


2.00 X (T/300)^°-4° 


2 


120) 


H2CO + + 


e- CO + H + H .... 


5.00 X (T/300)-°-50 


2 


121) 


HCO2 + + 


e- -> CO2 + H 


0.60 X (T/300)-°-S* 


2 


122) 


HCO2 + + 


e- -5> CO + + H .... 


0.84 X (T/300)-°-S4 


2 


123) 


H3CO + + 


e- H2CO + H 


2.00 X (T/300)^°-50 


1 


124) 


H3CO + + 


e-^CO + H + H + H 


2.00 X (T/300) 


1 


125) 


CO2 + + e 


- ^ CO + 


3.80 X (T/300)-''-50 


2 


126) 


CN + + e- 


-)-C-|-N 


1.80 X (r/300)-° '50 


2 


127) 


H2NCO + + e- HNCO + H .... 


5.0 X (r/300)-°-50 


1 


128) 


H2NCO + + e- -)■ H2 + NCO 


5.0 X (r/300)-0-60 


1 
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Table A7. Neutral-Neutral Reactions 



no. 




Chemical rea<;tion 


rate coefRcient(10 ^^cm^s ^) 


Reference 


129) 





+ OH -5> H + O2 


3.5 


2 


130) 


O 


+ CH ^> CO + H ... 


6.6 


2 


131) 


O 


+ CH2 —>■ OH + CH . 


49.8 X exp(— 6000/r) 


2 


132) 


O 


1 "IVTTT AT/^ 1 TT 

+ NH — !> NO + H ... 


11.6 


2 


133) 


O 


+ NH2 OH + NH 


1.39 X exp (-40/T) 


2 


134) 


o 


1 TVTTT TT"M/^ 1 TT 

+ NH2 — > HNO + H 


4.56 X exp(10/T) 


2 


135) 


o 


+ CN ->• CO + N ... 


4.36 X (T/300)° '*® exp (-364/T) 


2 


136) 


o 


+ CH3 H2CO + H 


13.0 


2 


137) 


o 


+ HNO ^ OH + NO 


6.0 


2 


138) 


o 


+ NCO NO + CO 


9.43 X (T/300)-0 09exp(-100/r) 


2 


139) 


c 


+ OH CO + H ... 


10 


2 


140) 


c 


+ NO ^ CO + N ... 


9 


2 


141) 


c 


+ NH ^ CN + H . . . . 


12 


2 


142) 


c 


+ NCO ^ CN + CO 


10 


2 


143) 


N 


+ OH ^ NO + H ... 


7.5 X (T/300)-° i8 


2 


144) 


N 


+ CH CN + H .... 


16.6 X (T/300)-° °8 


2 


145) 


N 


+ NH N2 + H .... 


4.98 


2 


146) 


N 


+ NO ^ N2 + .... 


3.75 X exp(-26/T) 


2 


147) 


N 


+ NCO ^ N2 + CO . 


4.0 


1 


148) 


H 


+ NCO ^ NH + CO 


12.6 X exp (-515/r) 


2 


149) 


CN + O2 -> NCO + . 


1.86 X (T/300)-° " X exp(40/T) 


2 
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Table A8. Reactions excited by shock heating part 1 



no. 




Chemical reaction 


rate coefRcient(cm^s ^) 


Reference 


150) 


Ha 


+ O^OH + H 


3.14 X 10-"(T/300)2'7 cxp (-3150/T) 


2 


151) 


Ha 


+ OH -5> HaO + H 


2.05 X 10-i2(T/300)i -'^2exp(-1736/T) 


1 


152) 


Ha 


+ HaO -5> OH + H + Ha 


5.8 X 10-9exp(-52900/T) 


2 


153) 


Ha 


+ C^CH + H 


6.64 X 10-1° exp (-11700/T) 


2 


154) 


Ha 


+ CH CHa + H 


5.46 X 10-10 exp (-1943/T) 


2 


155) 


Ha 


+ CHa ^ CH3 + H 


5.18 X 10-"(T/300)° i'^exp(-6400/T) 


2 


156) 


Ha 


+ CH3 ^ CH4 + H 


6.86 X 10-"(T/300)2 'i'4exp(-4740/T) 


2 


157) 


Ha 


+ N-5>NH + H 


1.69 X 10"-^cxp(-18095/T) 


2 


158) 


Ha 


+ NH -5> NHa + H 


5.96 X 10-" cxp (-7782/T) 


2 


159) 


Ha 


+ NHa NHs + H 


2.05 X 10-i'"''(T/300)^ *3cxp(-1400/T) 


2 


160) 


Ha 


+ NH3 ^ NHa + H + Ha 


1.5 X 10-** cxp (-42400/T) 


1 


161) 


Ha 


+ CN HON + H 


4.04 X 10-i=^(T/300)2''^'^cxp(-820/r) 


2 


162) 


Ha 


+ NH3+ NH4+ + H . . 


3.36 X 10-1'' exp (35. 7/r) 


2 


163) 


H + OH O + Ha 


6.99 X 10-"(T/300)2-80exp{-1950/T) 


2 


164) 


H + HaO OH + Ha 


1.59 X 10-ii(T/300)i-20exp(-9610/T) 


2 


165) 


H - 


1- CH ^ C + Ha 


1.31 X 10-1° exp (-80/T) 


2 


166) 


H - 


1- CHa ^ CH + Ha 


6.64 X 10-11 


2 


167) 


H - 


1- CH3 ^ CHa + Ha 


1.00 X 10-1° exp (-7600/T) 


2 


168) 


H - 


F CH4 ^ Ha + CHs 


5.94 X 10-"{T/300)-'-"° cxp (-4045/T) 


2 


169) 


H - 


1- NH ^ N + Ha 


1.73 X 10-ii{T/300)°-'''"cxp(-2400/T) 


2 


170) 


H - 


1- NHa -> NH + Ha 


5.25 X 10-i2(T/300)°-'^'^cxp(-2200/T) 


2 


171) 


H - 


1- NH3 NHa + Ha 


7.8 X 10-i3(T/300)2 -i0exp(-4990/T) 


2 
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Table A9. Reactions excited by shock heating part 2 



no. 




Chemical reaction 


rate coefRcient(cm^s ^) 


Reference 


172) 


H + 


H2CO -> HCO + H2 


4.85 X 10-12(7-/300)1-90 cxp(-1379/T) 


2 


173) 


H + 


HNO -5> NO + H2 


4.5 X 10-11 (T/300)"-'^2(jxp(-329/T) 


2 


174) 


H + 


O2 -5> + OH 


2.61 X 10-i"cxp(-8156/T) 


2 


175) 


O2 H 


h N ^ NO + 


2.26 X 10-i2(T/300)°-**'exp(-3134/T) 


2 


176) 


O2 + NCO CO + NO2 


8.1 X 10-11 exp(-773/T) 


1 


177) 


+ N2 -> N + NO 


2.51 X 10-i°exp(-38602/T) 


2 


178) 


+ 


H2O ^ OH + OH 


1.85 X 10-ii(T/300)''-95exp(-8571/T) 


2 


179) 


+ 


HON -5> OH + CN 


6.21 X 10-l"cxp(-12439/T) 


2 


180) 


+ 


NO ^ N + O2 


1.18 X 10-11 cxp (-20413/T) 


2 


181) 


+ 


CH4 ^ OH + CH3 


2.29 X 10-i2(T/300)2-2t' cxp (-3820/T) 


2 


182) 


+ 


NH3 -> OH + NH2 


1.89 X 10-11 cxp (-4003/r) 


2 


183) 


N + 


CN -5> N2 + C 


3.00 X 10-1'' 


2 


184) 


N + 


CO2 ^ NO + CO 


3.20 X 10-13 exp(-1710/T) 


2 


185) 


OH + CO CO2 + H 


2.81 X 10-i3exp(-176/T) 


2 


186) 


OH 


+ OH H2O + H 


1.65 X 10-i2(T/300)i i''exp(-50/T) 


2 


187) 


OH 


■f NH3 NH2 + H2O 


1.47 X 10-i3(T/300)2 05exp(-7/T) 


2 


188) 


OH 


■f NH2 NHa + 


2.26 X 10-l2(T/300)'' ''*'exp(-262/r) 


2 


189) 


OH 


+ CH4 -5> H2O + CH3 


3.77 X 10-"(T/300)2-42 exp (-1162/T) 


2 


190) 


OH 


+ NO ^ NO2 + H 


5.2 X 10-12 cxp (-15100/r) 


2 


191) 


Hc+ 


+ CH4 CH4+ + He 


5.10 X 10-11 


2 


192) 


Hc+ 


+ CH4 ^ CH2+ + He + H2 ... 


9.5 X 10-1° 


2 


193) 


He+ 


+ CH4 ^ CH3+ + He + H .... 


8.5 X 10-11 


2 


194) 


He+ 


+ CH4 ^ CH+ + He + H2 + H 


2.4 X 10-1° 


2 


195) 


He+ 


+ CH4 CHa + He + H+ .... 


4.8 X 10-1° 


2 


196) 


He+ + HNO -> NO + He + H+ 


1.0 X 10-s X (T/300)-°-5 


2 


197) 


He+ 


+ HNO NO+ + He + H 


1.0 X 10-8 X (T/300)-0-5 


2 



Table AlO. Reactions catalyzed by dust grains 



no. 


Chemical reaction 


rate coefficient (cm^s i) 


Reference 


198) 


H + H H2 


2.5 xlQ-i'' 


1 
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Table All. Average relative concentration data of species 



species 


mean relative concentration 


standard deviation of relative concentration 


initial relative concentration 


N+ 


8.530569 X 10"!° 


6.529162 X 10-° 


3.729334 x IQ-n 


He 


1.398994 X 10-1 


4.03005 X 10-5 


1.4 X 10-1 


NH+ 


1.479286 X 10"" 


1.5110507 X 10-11 


1.42177 X 10-11 


Na+ 


6.771889 X 10"* 


5.185934 X 10-° 


7.872506 X 10-* 


NO+ 


4.384701 X 10-12 


5.8279299 x IQ-i^ 


6.514852 X 10-12 


CUj 


2.()712():5 X 10 ^- 


1.51815:5 X 10 " 


1.07()720 X 10 ^- 


NH+ 


6.958145 X 10"" 


7.035392 X 10-11 


6.595174 X 10-11 


HON 


1.282398 X 10"* 


2.189310 X 10-° 


1.127323 X 10-* 


CH 


3.089765 X 10" 1° 


5.329314 X 10-1° 


2.247113 X 10-1° 


NCO 


7.930698 X lO"" 


7.993701 X 10-° 


8.471728 X 10-° 


CO2 


9.314529 X 10^1" 


7.334628 x 10^1° 


6.407260 X 10-1° 


N2 


1.017955 X 10-5 


8.040188 X 10-'^ 


1.206899 X 10-5 


CH+ 


1.812408 X 10^" 


1.762797 X 10-11 


2.137332 X 10-11 


H2O 


6.833447 x lO"'^ 


2.910295 X 10-'^ 


4.610721 X 10-'^ 


H3O+ 


6.044944 X 10"^ 


6.726546 x 10- 1° 


5.794733 x 10-i° 


NH3 


5.348606 X 10"* 


1.581378 X 10-8 


4.958232 x 10-* 


NH2 


1.794842 X 10"* 


1.925119 X 10-* 


1.773805 X 10-° 


11 


:-!.<s82()0S X 10 


:-!. 187012 X 10 


5,:-! 10500 X 10 


OH 


2.75309 X 10-8 


5.30112 X 10-* 


2.047385 X 10-* 


NO 


1.519176 X 10"'^ 


1.552755 X 10-'^ 


2.267862 x lO-'' 


N 


1.99560 X 10-« 


4.254789 x IQ-"^ 


3.036756 X 10'^ 


HCO+ 


8.217715 X 10-1° 


7.474217 X 10-1° 


8,544289 x 10-1° 


CO 


6.038563 X 10^^ 


9.499573 x IQ-^ 


6.27347 x 10-5 


H2 


4.980829 X lO^i 


3.248545 x IQ-^ 


4.980191 x 10^1 


Na 


3.42078 X 10-8 


5.17178 X 10-° 


2.327494 x 10-* 


ON 


7.377714 X 10-11 


6.538087 X 10-" 


1.708823 X 10-11 


C 


9.874322 x 10-* 


2.569366 x 10-* 


1.481840 X lO-'' 


H+ 


2.651750 X 10-1° 


4.785956 x 10-1° 


2.636774 x 10-1° 




1.160548 X 10-9 


1.658192 X 10-° 


9.953887 x 10-1° 


He+ 


5.963769 x 10-1° 


1.021091 X 10-° 


4.988382 x 10-1° 


O2 


2.297827 x IQ-^ 


3.048963 X IQ-^ 


1.793247 X 10-5 


C+ 


4.391799 X 10-9 


6.69398 X 10-° 


4.871363 X 10-° 





6.863861 X 10-5 


7.060412 X 10-'^ 


7.667686 x 10-5 


e~ 


7.665264 x 10-* 


1.741773 X 10-* 


8.70602 X 10-* 
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APPENDIX B: 
ESTIMATES 



MICROSCOPIC PHYSICS 



By using the Lagrangian particles, we can cleanly separate 
the macromixing due to turbulence from the micromixing 
due to the molecular diffusivity, an enormous advantage to 
this scheme, in contrast to the Eulerian schemes, in which 
these two effects are necessarily intermingled. The section 



9.8 of review paper by Veynante & Vervisch (20021 crys 



tallizes the distinction we have made in micromixing and 
macromixing. In the present Lagrangian formulation we 
have neglected the micromixing at the lower scales in the 
flow may be modelled as subsonic and incompressible. We 
also derive simple estimates for the microphysical transport 
coefficients. 

Next, we estimate the magnitudes of the mass diffusiv- 
ity coefficient D, and the specific turbulent kinetic energy 
dissipation rate e, in order to determine the critical turbu- 
lent mixing time (D/e)^^^. The neutral-neutral diffusivity 
Dun is given by the product of sound speed and the mean- 
free path length of neutral-neutral collisions. We denote r-Hj 
to be the intermolecular distance between two hydrogen 
molecules. We denote r^j^^ is the finite distance at intermolec- 
ular potential is zero. From Lennard- Jones theory, we take 
r]^^ — 2.5ri^^ to be the cut-off distance beyond which inter- 
molecular forces reduce to zero ( |Haile|1992| . Thus two neu- 
tral hydrogen molecules are attracted for r^j^^ < < ilii^- 
The collisional cross-sectional area for neutral-neutral colli- 
sions becomes (Tnn = tt (vii^ — V112) • Taking rj^^ — 2.93 A 
i Hirschfelder et al.|l964| ), ann = 6.07 x 10~^^cm^ 

The fiducial value of neutral-neutral diffusivity is given 

by 



(Bl) 



2.02x10" 



0.244 km s-i 



L 
1 pc 



ergg 



-\(B5) 



While all our estimates here will assume e is independent of 
length scale, turbulent intermittency results in fluctuations 
of this value ( |Pope|2000 1. 



1.00 x 10 



0.244 X 10* km s 



-)( 



4 X 10* cm-3 
1 



.6.07 X 10-i5cm2 

where rin is the neutral number density, and Unn is the cross 
section for neutral- neutral H2 molecular collisions. 

The neutral-neutral dynamic viscosity Unn is 

= PnCnAnn/2 



'(B2) 



Mnn ~ PnCn/ {2nnO'nn) • 



(B3) 



where pn is the neutral mass density, Ann is the mean-free 
path length for neutral-neutral collisions, and Cn is the sound 
speed for neutral species. The factor of 1/2 appears from 
angle-averaging over all distances between and Ann, and 
that their average particle velocities change linearly with dis- 
tance, when Ann is small. For a completely neutral medium, 
the kinematic viscosity is given by 



A^nn 
pn 



5.00 X 10' 



.0.244 km s- 



4 X 10* cm-3 
1 



6.07 X 10-^^ cm2 



-\(B4) 



These estimates have been used in calculating the Reynolds 
number in section [2^ 



The mean speciflc turbulent kinetic energy dissipation 
rate is 
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APPENDIX C: ANALYSIS OF DAMKOHLER 
NUMBERS OF REPRESENTATIVE SPECIES 

In this appendix, we identify the dominant chemical reac- 
tions for the Damkohler number calculations given in section 

The chemical network used in our codes have been tabu- 
lated in Appendix|A] The network consists of both one-body 
reactions, such as ionization by cosmic ray protons, and two- 
body reactions, including radiative recombination, charge 
transfer, recombination and neutral-neutral reactions. We 
have used the Damkohler number to classify the chemical 
evolution in the post-shock evolution of various species, in- 
cluding H2O, CH2, and CO in section ijs] We now derive ex- 
pressions for these Damkohler numbers which we employed 
in the text. 

To begin with, we consider H2O. The dominant creation 
reaction for H2O is 



The rate of change of concentration of methylene radical is 
then given by 



d[CH2] 
dt 



= fcii8[CH+][e-]-fci66[H][CH2] 



(C8) 



The Damkohler number for CH2 is then given by 



Da(CH2) = 



cooling 



^chemical 

((1.5[He] + 2.5P2I) X kB(Tp - Tan.biont)) 



A(n,r) 
fcii8[CH+][e- 



fcl66[H][CH2 



[CH2 



(C9) 



110 > HgO^ 



0.295 H2O + 0.705 OH + 1.705 H ; ku 



The numbering of reactions is the same as that in Appendix 
[A} and k denotes the rate constant of a reaction, with the 
subscript denoting the reaction number in Appendix]^ This 
convention is carried forward throughout the Appendix [C] 
Similarly, the dominant destruction reactions for H2O are 



42 > C+ + H2O HCO+ + H 
50 > HCO+ + H2O ^ H3O+ + CO 
61 > Ut + H2O ^ H3O+ + H2 



^42 



We can insert typical values of physical variables in the equa- 
tion [C9] to get an expression of the scaling of the Damkohler 
number of CH2 under a representative shock. We select a 
shock (a part of the trajectory is shown in figure [2]| with 
jVIach number of 3.5. The post-shock temperature Tp is then 
calculated to be 27.9 K. The actual concentration of each 
species is indicated in the denominator in the scaled equa- 
tions which follow below. 

To connect the Damkohler number to the post-shock 
temperature, we take the ambient temperature to be 10 K. 
Then the post-shock temperature Tp is given by the 
Rankine-Hugoniot jump condition 

Tp _ {27M^-(7-l)}{(7-l)A1^+2} ^^^^^ 



(C2) 



This result is used in figure [6] and in section !j5] 

The rate of change of water at any given instant, taking into 

account the dominant reactions, is 



= 0.295fciio[H30+][e-] - fc42[C+][H20] 
-A;5o[HCO+l[H2 0] - fceip+lpaO] 



(C3) 



Using the above information, the chemical time scale for 
water to get to equilibrium is given by 

[H2O] 



tchcm (H2O) 



d[H2 0] /dt 



(C4) 



The Damkohler number for water Da(H2 0) is hence given 
by (using the definition in equations |45[ [C4] and section [sTs] ) 

Da(H20) = ^'^""""^ 

tch cmical 



(C5) 



Next, we consider methylene radical, which we pre- 
sented as an example of rapidly varying (shock-enhanced) 
species in ^ We now derive the expression used for calculat- 
ing its Damkohler number. The dominant creation reaction 
for methylene radical CH2 is 



118 > CHJ + e~ CH2 + H ; fcn 



(C6) 



while the dominant destruction reaction for methylene rad- 
ical CH2 is 



166 > H + CH2 ^ CH -f H2 : fci 



(C7) 



Tambi' 



(7 + 1) A1 = 



Using the above relation, we can predict the range of the 
Damkohler number for methylene radical as a function of 
time, provided we know the pre-shock concentration. Signif- 
icantly, this also implies that the Damkohler number is not 
an universal property, and can vary if the chemical concen- 
tration changes. In the limit of a strong shock we obtain 



Tp _27(7-l)X^ 



(7+1)' 



(Cll) 



On assuming a completely molecular gas with 7 = 7/5, we 
obtain 



To 



^ambic 



Ta, 



7M^ 
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We have presented carbon monoxide as an example of 
chemically frozen species in !|5] The dominant creation re- 
action for carbon monoxide CO is 



108 > HCO+ -f e" ^ CO + H ; fci 



(C13) 



while the dominant destruction reactions for carbon monox- 
ide CO are 



21 > He+ + CO ^ C+ + O + He ; kii 



57 > Ht + CO ^ HCO^ 



H2 



(C14) 



Based on the above reactions, the rate of change of carbon 
monoxide is given by 



d[CO] 
dt 



= fcios [HCO+] [e-]-fc2i [CO] [H+]-fc57[CO] [He+] .(C15) 



Using the above relations, the Damkohler number for carbon 
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monoxide is given by 

^cooling 



Da(CO) 



tchemical 

((1.5[He] + 2.5[H2]) x kB(Tp - T^^uent)) 
A(n,T) 

, fci08[HCO+][e-] - fe2i[C0][H+] - fc57[CO][He+] ^^^^^ 



